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As is well-known, in the conventional formulation of Bogoliubov's theory of an interacting Bose 
gas, the Hamiltonian H is written as a decoupled sum of contributions from different momenta of 
the form H = J^k^o where each Hamiltonian Hk describes the interaction of bosons in the 
condensed k = state with bosons in the momentum modes ±k. Then, each of the single-mode 
Hamiltonians Hk is diagonalized separately, and the resulting ground state wavefunction of the 
total Hamiltonian H is written as a simple product of the ground state wavefunctions of each of the 
single-mode Hamiltonians Hk- We argue that, while this way of diagonalizing the total Hamiltonian 
H may seem to be valid from the perspective of the standard, number non-conserving Bogoliubov's 
method, where the k = state is removed from the Hilbert space and hence the individual Hilbert 
spaces where the Hamiltonians {-Hk} are diagonalized are disjoint with one another, from a number- 
conserving perspective this diagonalization method may not be adequate since the true Hilbert 
spaces where the Hamiltonians {-ffk} should be diagonalized all have the k = state in common, 
and hence the ground state wavefunction of the Hamiltonian H may not be written as a simple 
product of the ground state wavefunctions of the HkS. In this paper, we give a thorough review of 
Bogoliubov's method, and discuss a variational and number- conserving formulation of this theory 
in which the k = state is restored to the Hilbert space of the interacting gas, and where, instead of 
diagonalizing the Hamiltonians Hk separately, we diagonalize the total Hamiltonian J? as a whole. 
When this is done, the spectrum of excitations of the system changes from a gapless one, as predicted 
by the standard, number non-conserving Bogoliubov theory, to one which exhibits a finite gap in 
the k — > limit. 



I. INTRODUCTION 



There has been a lot of interest in the properties of 
interacting Bose systems during the past fifteen years 
following the experimental observation of Bose-Einstein 
condensation (BEC) in ultracold vapors of trapped 
rubidium- and sodium gasesi^ On the theoretical front, 
much of the effort to understand the properties of these 
systems has focused on using approaches of the Bogoli- 
ubov type at the lowest temperatures ^ — and finite tem- 
perature generalizations of Bogoliubov's theory^— at 
higher temperatures. Yet, and despite considerable ac- 
tivity, a closer look reveals that several aspects of Bo- 
goliubov's theory are still not fully understood^r— In 
this paper, we will examine Bogoliubov's theory in quite 
some detail, and raise a few issues which have been over- 
looked in the past, and which in the author's view still 
need to be clarified. Chief among the questions that we 
want to address is the question of trying to understand 
how the non-conservation of particle number, which is 
a landmark feature of Bogoliubov's theory, and the de- 
coupled way in which the Hamiltonian is diagonalized, 
affect the nature of the Bogoliubov ground state. In- 
deed, and as is well-known, in the standard formulation 
of Bogoliubov's theory, the creation and annihilation op- 
erators of the condensate, and ao respectively, are re- 
placed by the c-number ~ vOV, N being the total num- 
ber of bosons in the system. The Hamiltonian H is then 
written as a decoupled sum of contributions from differ- 
ent momenta of the form H = X^k^o where each 
Hamiltonian H k describes the interaction of bosons in 
the condensed k = state with bosons in the momen- 



tum modes ±k. Then, each of the single-mode Hamilto- 
nians i?k is diagonalized separately and the ground state 
(GS) wavefunction of H is obtained as the product of 
the GS wavefunctions of the H^s. In this paper, we 
shall argue that, while this way of diagonalizing the to- 
tal Hamiltonian H may seem to be valid from the per- 
spective of the conventional, number non-conserving Bo- 
goliubov's method, where the k = state is removed 
from the Hilbert space and hence the individual Hilbert 
spaces where the Hamiltonians {H k } are diagonalized 
are disjoint with one another, from a number-conserving 
perspective this diagonalization method is not appropri- 
ate since the true Hilbert spaces where the Hamiltonians 
{Hk} should be diagonalized all have the k = state 
in common, and hence the ground state wavefunction of 
the Hamiltonian H may not be written as a simple prod- 
uct of the ground state wavefunctions of the H^s. We 
then shall discuss a variational, number-conserving gen- 
eralization of Bogoliubov's theory in which the k = 
state is restored into the Hilbert space of the interacting 
gas, and where, instead of diagonalizing the Hamiltoni- 
ans _ffk separately, we diagonalize the total Hamiltonian 
H as a whole. When this is done, the spectrum of ex- 
citations of the system changes from a gapless one, as 
predicted by the standard, number non-conserving Bo- 
goliubov method, to one which exhibits a finite gap in 
the k — > limit. 

The rest of this article consists of two main parts. The 
first part is almost entirely devoted to gaining a deeper 
understanding of Bogoliubov's theory, and consists of 
Sees. |H] through [VlJ In the second part of the paper, 
consisting mainly of Sec. I VIII we go beyond the conven- 
tional formulation of Bogoliubov's theory by trying to en- 
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force the conservation of the number of bosons in the sys- 
tem, making sure we keep an accurate count of the num- 
ber of bosons in the k = state. As mentioned above, 
when this is done, the results of Bogoliubov's theory are 
changed both in quantitative and qualitative ways, and 
these changes are discussed in detail in Sec. I Villi 

We now want to give the reader a quick overview of 
the contents of each one of the remaining eight Sections 
of this manuscript. In Sec. m we shall start by review- 
ing the standard Bogoliubov treatment of an interacting 
Bose gas. As is well-known, in this treatment the total 
number of bosons is not conserved, which leads to a num- 
ber of unphysical features. As an example, the average 
number of bosons Nk in the single-particle state of mo- 
mentum k is found to diverge like 1/k as k — > 0, which 
does not of course make much sense for a system with a 
fixed number of bosons N. In Sec. Mil we discuss this 
and other similar unphysical predictions of Bogoliubov's 
method, and their consequence on the evaluation of ex- 
pectaion values of physical observables in the Bogoliubov 
ground state. To be able to formulate Bogoliubov's the- 
ory within a number-conserving framework and so get rid 
of the aforementioned unphysical features, in Sec. HVlwc 
proceed to derive a number-conserving version of Bogoli- 
ubov's Hamiltonian, which we show can be written as a 
sum of decoupled Hamiltonians for each momentum 
mode k, H = X^k^o^k- Then, in order to substantiate 
the claim made above according to which each one of the 
single-mode Hamiltonians corresponding to a given 
momentum mode k is diagonalized independently in Bo- 
goliubov's theory, in Sec. |V]we present a detailed anal- 
ysis of a variational, number-conserving approach to the 
single- mode Hamiltonian H^, and show explicitly that 
the results obtained from the analysis of this single-mode 
Hamiltonian perfectly coincide with the results of Bogoli- 
ubov's method in the thermodynamic N — > oo limit. At 
this point, we will find it useful to try to assess the ac- 
curacy of the above variational analysis, which is why in 
Sec. IVII we perform an exact numerical diagonalization 
of the single-mode Hamiltonian studied in the pre- 
ceding Section. Our numerical results corroborate the 
variational approach of Sec. [Vj which is found to be sur- 
prizingly accurate. The main conclusions of the analysis 
carried out in Sees. [V] and IVII is that the treatment of 
the boson Hamiltonian carried out in Bogoliubov's theory 
amounts to diagonalizing each one of the Hamiltonians 
Hk (representing the kinetic energy plus the interaction 
energy of bosons having a momentum +k or — k with 
the condensate) independently from the other momen- 
tum contributions -ffk'(^k)- While this would be perfectly 
legitimate if the various Hubert spaces used to diagonal- 
ize the Hamiltonians were disjoint, in our case all 
these Hilbert spaces have the k = state in common, 
and so the above decoupled diagonalization procedure, 
strictly speaking, is mathematically not valid. 

In the second part of the paper, and in order to assess 
the quantitative accuracy of diagonalzing the -ffk's in the 
decoupled way described above, we generalize the varia- 



tional approach of Sec. [V] for the single-mode Hamilto- 
nian i?k to the total Hamiltonian H = X)k^o ^ k m ^ ec - 
IVIII and we do that in such a way as to keep an accurate 
count of the number of bosons in the k = state, and 
with the requirement that the total number of bosons 
in all momentum modes be conserved. This leads, quite 
surprizingly, to an excitation spectrum of bosons which 
presents a gap as k — > 0, by contrast to the usual Bo- 
goliubov method in which the spectrum of excitations 
is gapless. In Section IVIIII we give a discussion of our 
results, and in Section HXl we present our conclusions. 

II. REVIEW OF THE STANDARD, NUMBER 
NON-CONSERVING FORMULATION OF 
BOGOLIUBOV'S THEORY 

We shall begin by reviewing the standard, number non- 
conserving formulation of Bogoliubov's theory. To this 
end, let us consider the Hamiltonian of a gas of N spinless 
bosons, interacting through a two-body potential v(r): 

H = J dr*t(r)(-— J*(r) 

+ i J dvdr' *t(r)* + (r')w(r - r')*(r')*(r), (1) 

where ^(r) is a boson field operator in second quan- 
tized language, which has the Fourier decomposition 
^(r) = ^ k ake lk r / VV (periodic boundary conditions 
will be assumed throughout this paper), ak being a bo- 
son annihilation operator, and V being the volume of 
the system. In Fourier space, the Hamiltonian H has the 
following expression: 

H = e k a k a k + ^7 J2 J2 K<l)4+q4'-q a k'ak, (2) 

k^O k : k' q 

where w(k) is the Fourier transform of v(r) and £k = 
H 2 k 2 /2m is the kinetic energy of bosons with wavevector 
k. 

According to Bogoliubov's prescription^^ one re- 
places the creation and annihilation operators aj and do 
of the k = state by v^Oj where Nq is the number of 
bosons in the condensate, upon which the above Hamil- 
tonian takes the form H = H Q + H 2 + H3 + H 4 , with 
(n = No/V is the density of condensed bosons): 

Ho = \vn 2 v(0), (3a) 

H 2 = ^ |Cka k fl k + ^ov(k) [a k a -k + a ka-k] }, (3b) 

k^O 

^ 3 = V V ^ V ^ [ a k+q ak0t l + a i< a q a k+q] . (3c) 
q,k 

^ = ij E w(q)4 +q 4,_ q a k /ak. (3d) 

q,k,k' 
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In Eq. (13b[) . £ k is the shifted boson energy: 



£k 



•no[«(k)+«(0)], 



(4) 



and it is understood that no creation and annihilation 
operators of the condensate (aj and a ) appear on the 
rhs of Eqs. (3c]) and (3d)) . 

As it can be seen, the new expression of the Hamilto- 
nian after the Bogoliubov prescription is performed does 
not conserve the number of bosons. To deal with this un- 
physical artifact of Bogoliubov's prescription, Bogoliubov 
suggests that one should work with the "grand canonical 
Hamiltonian" H = H — /j,N, \i being the chemical poten- 
tial and the operator N — J^ k a k a k being the total num- 
ber of bosons. This amounts to replacing Hq and Hi in 
Eqs. by H = Ho — M^o and H 2 = H 2 — (iNi, respec- 
tively, where Ni = X^k^o a ik ak * s ^ ne operator counting 
the total number of bosons outside the condensate. The 
ground state energy is then found by diagonalizing H2, 
which is done by introducing a new set of creation and 
annihilation operators, a k and ak, such thai 



26. 



a k = u k a k + fka_ k , 



a k = u k a k 



v k a_ k . 



(5) 



These expressions can easily be inverted, with the result: 
a k = u k a k - i> k a^ k , a k = M k a k - t> k a- k . (6) 

In Eqs. (5J and ©, u k and v-^ are assumed to be real, 
spherically symmetric functions of the wavevector k. For 
the newly defined operators a k and a k to describe boson 
excitations, it is required that they should obey bosonic 
commutation relations: 

[ak, c*k'] = [a k , a k ,] = 0, [a k , a£,] = £ k ,k', (7) 

which results in the condition 26 



22-1 
tlv — = 1. 



(8) 



If we now use Eq. (5]) to rewrite the Hamiltonian H2 in 
terms of the operators a k and a k , and if in addition we 
require that the resulting expression not contain products 
of the form a k a_ k or a k al k , we obtain the following 
expressions of the so-called "coherence factors" w k and 
w k : 



1 /6c - M 



( 



1 



1 /£k - M 



Ek / 2 V _E k 

where _E k is the energy spectrum: 



i?k-A/(ek-M) 2 -^(k)2. 



-1 , (9) 



(10) 



The Hamiltonian H2 then takes the quadratic form: 

H 2 = E 2 + Y, E k ala k7 (11) 

k^O 



where we denote by E2 the following quantity: 

^ 2 = \Y1 [E k -e k -n v(k)]. (12) 

k^O 

Requiring the spectrum £" k of Eq. (fTOf to be gapless^ as 
k — > results in the following expression of the chemical 
potential: 

fi = n v(0), (13) 
upon which E^ takes the celebrated Bogoliubov form: 



£ k = yek[ek + 2no«(k)]. (14) 
Using Eq. (fl"3"j) into Eq. ©, one finds: 

2 l/e k + n w(k) \ 2 l/e k + n v(k) 
Uk= 2l EL + 1 > Uk= 2l El 1 



(15) 

The Bogoliubov ground state (BGS) of the interacting 
bose gas, which we shall denote by the symbol \^b), is 
defined as the vacuum state for the a k operators, and 
satisfies the condition: 



a k |W B )=0, forallk^O. 



(16) 



With this definition, it is easy to verify that E 2 is the 
expectation value of the Hamiltonian H2 in the BGS, i.e. 
E 2 = (* b \H2\*b)- 

We now pause for a moment to emphasize the distinc- 
tion which is usually done in Bogoliubov's theory between 
the operators a k and a k on one hand, and the operators 
a k and a k on the other. For while the former create and 
annihilate actual bosons, the latter are thought to cre- 
ate and annihilate "quasiparticles" , which are identified 
as collective "phonon" modes, mainly because the corre- 
sponding excitation spectrum _E k has a linear behavior 
at long wavelengths: 



shk, 



where 



= y/n B v(0)/m 



(17) 



(18) 



is identified as the speed of sound (here and in the rest 
of this paper, tib = N/V is the density of bosons). We 
will come back to this identification in Sec. IVIII Cl below. 
For the moment, we want to have a more detailed look at 
a number of problematic features of Bogoliubov's theory. 
These being more or less known features of this method, 
our main objective here is to provide a motivation for 
the variational formulations of Bogoliubov's theory that 
will be studied in Sees. IVl through IVIII This will be the 
subject of the following Section. 
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III. ISSUES WITH BOGOLIUBOV'S THEORY 

As mentioned above, in this Section we want to discuss 
a few questionable aspects of the number non-conserving 
formulation of Bogoliubov's method as a way to motivate 
the variational treatment that will be presented in Sees. 
IVl and I VIII We will start by looking at the ground state 
expectation values of quadratic combinations of creation 
and annihilation operators. 

A. Divergence of the depletion Nk = ( a k a k) and of 
the anomalous average (aka-k) in the k — > limit 

The first quantity we will discuss is the expectation 
value: 

A^k = (4«k> = <*sl4 a k|*i3>, (19) 

which gives the average number of bosons in the single- 
particle state of momentum k. In the standard formula- 
tion of Bogoliubov's theory, this quantity is given by: 

N^ V l = \(^±^-l). (20) 

Using the expression of i? k given in Eq. (|14p , it is not dif- 
ficult to see that iVk diverges like 1/k as k — > 0, which is 
of course unphysical, given that we started from a system 
of fixed number N of bosons. 

The divergence of as k — > will have important 
ramifications, as it will affect the expectation value of 
one-body operators. Indeed, if we consider a given phys- 
ical observable which is described by a one-body opera- 
tor O of the form O = ^ k Ok.a^ak, then the expectation 
value of O in the Bogoliubov ground state is given by: 

(0)=J20kN k . (21) 

k 

Since the result (|20l) for Nk is unphysical near k = 0, we 
see that the result (|2 1 1) for the expectation value (O) in- 
cludes terms whose contribution is unphysical near k = 
as well. As an example, let us consider the Fock part of 
the Hamiltonian Hp = X)k^o n o u (k)a k a k- The contri- 
bution of any particular k mode in this last expression 
to the ground state energy is given by noi>(k)./Vk, which 
diverges like 1/k as k — > 0. This divergence cannot of 
course be taken at face value, because it originates in the 
unphysical divergence of the quantity Nk, which should 
be finite and bounded by the total number N of bosons 
in the system for all values of the wavevector k. 

A closely related problematic feature of Bogoliubov's 
theory has to do with the divergence of the anoma- 
lous averages (okfl-k) = s| a k a -k| v f' b) an< i ( a ik a -k) = 
(\E , s|a ] J i aL k |\E , s) in the k —> limit. Following an ar- 
gument made by de Gennes for the BCS ground state of 
superconductors^, one can argue that an average such as 



(ak<2-k) can be given a sense, from a number-conserving 
perspective, if it is understood as: 

(akO_ k ) - (*(JV - 2)|a k a_ k |*(7V)), (22) 

where we denote by |\&(iV)) the normalized ground state 
wavefunction of a system of N bosons. Written in this 
form, (a k a_k) can be interpreted^! as the probability 
amplitude for finding the system in the ground state 
\ty(N — 2)) when a pair of bosons (k, — k) is removed 
from the ground state ^(JV)). Following this line of 
thought, if we think of (aka_k) as a probability ampli- 
tude, then it should be a bounded quantity, and should 
not diverge for any value of the wavevector k. Unfortu- 
nately, in the standard Bogoliubov theory, the expression 
of the anomalous average is given by: 

/ t t n ov(k) 

(a k a-k) = ( a k a - k ) = ► ~°° as fc ^ 0, 

(23) 

and diverges (negatively) as k — > 0, which is of course 
incompatible with an interpretation of (a k a -k) ^ n terms 
of probabilities. The divergence of the anomalous aver- 
age (akO-k) will also have major consequences on the 
evaluation of expectation values of quantities involving 
quadratic products of operators of the form ak&_ k or 
a ik a -k- Again, as an example, if we consider the so- 
called "pairing" part of the Bogoliubov Hamiltonian, 
Hp = | J2k^o 7l oi'(k)(a k aLk + a kO-k), then the contri- 
bution of any particular term in Hp corresponding to a 
given wavevector k to the ground state energy is given by 
— [nof(k)] /2Ek, and again diverges like —1/k as k —> 0. 
This divergence is again not to be taken too seriously, 
since it originates solely from the unphysical divergence 
of the anomalous average (akfl-k) m Eq. ([23]) . We will 
look at these divergences in more detail below (see Sec. 
rVj) . For the moment, we want to have a closer look at 
the physical meaning of the operator H in Bogoliubov's 
method. 



B. Physical meaning of the operator H in 
Bogoliubov's method 

We now want to discuss the physical meaning of the 
Hamiltonian H in the standard formulation of Bogoli- 
ubov's method in view of the fact that the ground state 
energy of the system is the expectation value of the 
"grand-canonical" Hamiltonian H = H — fiN, and not of 
the Hamiltonian H itself. As an aside, we shall observe 
that the use of the term "grand-canonical Hamiltonian" 
to describe the operator H is quite misleading, given that 
the term "grand-canonical" has a very specific meaning 
in ensemble theory of statistical mechanics, and refers 
to a very specific statistical ensemble which only has a 
meaning at nonzero temperatures. We also would like 
to remind the reader that, in ensemble theory, various 
operators have the same physical meaning in the grand- 
canonical ensemble as in any other statistical ensemble. 
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In particular, the operator describing the total energy 
in the grand-canonical ensemble is the Hamiltonian H, 
not the combination H — fj,N, which in fact describes the 
free energy. It is therefore extremely surprizing that in 
Bogoliubov's theory the ground state energy is found by 
taking the expectation value of H = H — fiN instead of 
the expectation value of the Hamiltonian H itself^ 

In what follows, we want to go further and show that 
the use of the Hamiltonian H instead of H can even lead, 
in certain situations, to nonsensical results, which will al- 
low us to argue that the Hamiltonian H has no physical 
meaning in the standard formulation of Bogoliubov's the- 
ory. To this end, let us consider states that are generated 
by successive action of the "quasiparticle" operator a k on 
the BGS. It is easy to verify that these states are eigen- 
states of the "grand-canonical" Hamiltonian H 2 . For ex- 
ample, one can show that the normalized state defined 
by: 



I* 



(n)\ 



(24) 



and corresponding to the excitation of n "quasiparticles" 
of wavevector k is an eigenstate of Hq, with eigenvalue 
[E 2 +nE k ]: 



[E 2 +nE k ]\<S> 



Similarly, a state of the form: 



(25) 



(26) 



and corresponding to the excitation of n "quasiparticles" 
with distinct wavevectors ki,k2, • • • , k„, is an eigenstate 
of the grand canonical Hamiltonian H 2 with eigenvalue 
[E 2 + E kl + ■ ■ ■ + E kn ] : 



H 2 \$ kl ,..., kn ) = [E 2 +E kl 



ki,— .k,, 



(27) 



Let us focus for simplicity on the state with one "quasi- 
particle" of momentum k. This is the state given by: 



|*k>=ai|* fl ) 



(28) 



Such a state is believed to have an excitation energy 
which is gapless in the k — > limit. This belief originates 
in Eq. (|TH) . where E k is identified as the excitation en- 
ergy of the a k "quasiparticles" . Indeed, if we calculate 
the excess energy with respect to the Bogoliubov ground 
state of the state \Q k ) using the "grand canonical" Hamil- 
tonian H 2 , wc find: 

(* k |[flo + H 2 ]\$ k ) - (^ B \[H + H 2 ]\^ B ) = E k . (29) 

At this stage, it is useful to remember that the true 
quadratic part of the original Hamiltonian H is H 2 , not 
H 2 . Consequently, in evaluating the excitation energy 
(not the excitation free energy!) AE^ x ' a of a "quasipar- 
ticle" a k , one should use the Hamiltonian H 2 instead of 



H 2 , and define AE^' a as follows (we here use the super- 
script a to indicate that we are calculating the excitation 
spectrum of the a\. "quasiparticles" ) : 



($ k | [Ho +H 2 ] |$ k ) - | [H +H 2 ] | tt B ) . (30) 



Performing the calculation, we find: 



AE ex,a 



E k 
E k 



n u(0)[?i k + v ) 



^kJ' 



n v(0) 



e k + n v(k) 
y/ewisw + 2n v(k)} 



(31a) 
(31b) 



As it can be seen, the energy cost to excite a "quasipar- 
ticle" a k , when calculated using H 2 , is no longer given 
by the gapless expression E k . There is now an additional 
term, which furthermore diverges like 1/fc as k — > 0. We 
are therefore faced with the bizarre situation where the 
operator H 2 has no physical meaning of its own what- 
soever, and where the energy observable rather mysteri- 
ously is no longer represented by this last operator, but 
by the combination H 2 = H 2 — (aNi, which, in principle, 
represents the free energy of the system. Because of the 
non-conservation of particle number, the use of the true 
energy operator H 2 in the formulation of Bogoliubov's 
theory given in Sec. |H] can lead to non-sensical results, 
as we have seen above with the excitation spectrum of 
the a k "quasiparticles". 



C. Ad hoc nature of the canonical transformation 
from the a^'s to the a^'s 

We now turn our attention to another questionable as- 
pect of Bogoliubov's method, having to do with the phys- 
ical content of the canonical transformation from the ctk's 
to the ctk's, and the rather ad hoc fashion in which the 
commutation relations between the excitation operators 
a k and a k are introduced. Taking note of the fact that 
these operators, in the standard formulation of Bogoli- 
ubov's theory, do not conserve the number of bosons, be- 
low we will show that, in a number-conserving approach, 
the canonical commutation relations between the ak's 
are not at all exact, but only approximate, and are only 
valid if the depletion of the ground state is small. For sit- 
uations where this condition is violated, such as for liquid 
Helium at low temperatures, where the number of bosons 
in the condensate is only about 10% of the total number 
of bosons in the system, the canonical commutation re- 
lations are no longer valid, and should be replaced by 
more involved expressions. For more details, the reader 
is referred to Sec. [V] where the commutation relations 
between the a k s will be obtained from first principles, 
and not imposed a priori, within a number-conserving 
variational approximation. 
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D. Diagonalizing the various Hamiltonians 
corresponding to different momentum modes 
independently from one another 

Perhaps one of the less obvious criticisms one can make 
of Bogoliubov's method, and possibly one with the most 
far-ranging ramifications, has to do with the fact that, 
in this theory, the total "grand-canonical" Hamiltonian 
H = H — fiN is not diagonalized as a whole, but rather, 
each momentum contribution to H is diagonalized 
separately and independently from the other momentum 
contributions. One way to see this is by writing the to- 
tal Hamiltonian H as a sum of contributions from dif- 
ferent wavevectors k of the form H = J^k^o Hk- Then 
it is easy to convince oneself that the canonical trans- 
formation method of Sec. [TT] is in fact a diagonalization 
of each one of the Hamiltonians independently from 
the other Hamiltonians -ffk'(^k)- This process of diago- 
nalizing each of the single-mode Hamiltonians sepa- 
rately would be absolutely sound if the Hilbert spaces in 
which these Hamiltonians are diagonalized were totally 
disjoint from one another. This is unfortunately not the 
case here, with all these Hilbert spaces having the k = 
state in common. Of course, from the perspective of Bo- 
goliubov's method, where the dynamic properties of the 
condensate are ignored and the k = state is summar- 
ily removed from the Hilbert space, this is thought to be 
a perfectly legitimate way to proceed. However, in this 
paper, we shall challenge this view, and explicitly show 
that, when the conservation of particle number is prop- 
erly taken into account, and the k = is restored as an 
essential part of the Hilbert space used to describe the 
system, then the direct diagonalization of the full Hamil- 
tonian H gives very different results from diagonalizing 
each of the momentum components separately, as is 
done in the standard Bogoliubov formulation. 

Since the non-conservation of particle number in the 
standard formulation of Bogoliubov's theory is at the core 
of most of the conceptual difficulties discussed through- 
out this Section, it seems worthwhile to try to eliminate 
these difficulties by reformulating the theory within a 
number-conserving framework. This will be done in Sec- 
tions E] through ELD But, before we do so, we want 
to present a rather detailed derivation of Bogoliubov's 
Hamiltonian in a number conserving approach. This will 
be done next. 



IV. DERIVATION OF BOGOLIUBOV'S 
HAMILTONIAN WITHIN A 
NUMBER-CONSERVING APPROACH 

We now want to derive, within a number-conserving 
framework, a simplified version of the Hamiltonian ([5]) 
which is simple enough to allow for straightforward eval- 
uation of expectation values of physical observables, and 
at the same time captures the essential physics of a trans- 
lationally invariant system of interacting bosons at zero 



temperature. Throughout this paper, it will be assumed, 
without loss of generality, that the bosons are confined 
in a cubic box of size L, and we shall use periodic bound- 
ary conditions, implying that momenta will be quantized 
with the following wavevectors: 



2 k 



[tlx i Tly j Ylz ) 1 



(32) 



where the n^s are integers such that —oo<n x ,n y ,n z < 



oo. 



We shall start by specifying the Hilbert space we will 
use to describe our system. In the occupation number 
representation, a system of N bosons can in general be 
in any one of the states 



I rni-n 



= \N- 



(33) 



having ni bosons in the single-particle state of wavevec- 
tor kj and to, bosons in the single-particle state of 
wavevector — k; (for all i such that 1 < i < oo), and 
\N — X^i( n » + m i)] bosons in the single-particle state 
with wavevector k = 0. More formally, if we denote by 
|0) the vacuum state for bosons, then IV'™ 1 ---^ 00 ) is the 
ket defined by: 



K 1 .: 



/ t^-5:.=i(™-+ m ') 



n 



-|0>. 



(34) 



It then follows that the most general expression of the 
wavefunction |\I/(jV)) of a system of N bosons in the oc- 
cupation number representation is given by: 



E 



E 



s~im\ ,...,m 
^ n i ,...,noo 



I rrii"-n, 



(35) 



In the above equation, the C^yy^° 's are complex num- 
bers to be determined by diagonalization of the Hamilto- 
nian H of the system. The summations extend over the 
range < n^m, < N, subject to the constraint that the 
number of bosons in the state k = in any ket of the 
basis must be positive: 



./V — (ni + mi) ■ 



0>o. 



(36) 



Now, it can be easily verified that the total momentum 
operator P = ^ k Kka^a^ commutes with the Hamilto- 
nian H . This implies that H and P can be diagonalized 
simultaneously, and that eigenstates of H can be labeled 
by a definite value of the total momentum operator P. 
In particular, for a system of bosons at rest, the ground 
state is the translationally invariant state corresponding 
to the eigenvalue P = of the total momentum^ 



P|*(JV)) = 0. 



(37) 
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This condition imposes a rather cumbersome constraint 
on the coefficients C™ 1 ' '^™ 00 , namely: 

Tll,mi Moo .moo 

+ hk^inoo - moo)] IVCV-'.C ) = °- ( 38 ) 

The easiest way to satisfy the above constraint is to re- 
strict the Hilbert space to one where there are as many 
bosons with momentum ki as there are bosons with mo- 
mentum — ki, i.e. we choose to work in the restricted 
Hilbert space such that rii — rrii, for all values of the 
index i. Eq. (|35)l then becomes: 

MN)) = E ' ' ' E IVW- ,nJ (39) 

771 "cc 

where now iVW-.moo) is given by: 

iV'ni,- ,«oo) = l-^V - 2ni 2noo;ni,rai; . . . ; raoo^oo), 

(40a) 

= t a °) — 17 , ' I") ( 40b ) 

We now want to isolate those terms in the Hamiltonian 
which will give significant contributions to the ground 
state energy in states of the form (|39|) . For dcfinitcncss, 
let us rewrite the interaction part V of the many-body 
Hamiltonian, which is given by: 

^ = ^7 E E «(qK +q «l- q a k'ak. (41) 

k,k' q 

In the above expression, we shall first separate the 
Hartree terms, corresponding to k + q = k and k' q = 
k', i.e. to q = 0; and the Fock (or exchange) terms, 
which correspond to k + q = k' and k' q = k, i.e. to 
q = k' k. Hence, we shall write the interaction term V 
in the form: 

v = v H + % + — E E v (q) a k+ q 4<-q a k<ak, 

k,k' q#0,k'-k 

(42) 

where the Hartree and Fock contributions, Vh and Vf 
respectively, are given by: 

Vh = E 44 a k'Ok, (43a) 

k,k' 

= ^-N(N-l), (43b) 

^ = ^E E Kk' - k)4,4a k ,a k . (43c) 
k k'(#k) 

We next consider the so-called "pairing" terms. These 
are obtained by letting k' = — k both in Vf and in the 



last term of Eq. (|42|) . This allows us to write: 
V = V H + V F + V P 

+ ^yY, E v (q) a k +q a L-q a k'ak, (44) 

k k'(^-k) q#0,k'-k 

where we defined: 

^f=2vJ2 E v ( k ' -k)4'4ak/ak, (45a) 

k k'(^k,-k) 

= ^7 E E w (q)4 + q a -(k+q) a ka-k. (45b) 

k q#0 

It can be shown that the last term on the rhs of Eq. 
(|44| has zero expectation value in the state |\t(iV)) of 
Eq. (|39p. and hence it shall be discarded. Moreover, in 
keeping with the spirit of Bogoliubov's method, both in 
V F and Vp we shall only retain terms in which either k 
or k' is zero, an approximation which is believed to be 
valid if the depletion of the condensate is small. 29 Hence, 
we shall write: 

V = V H + V' F + Vp, (46) 

with: 

% ~ — E u(k)aja 4 a k, (47a) 

k/O 

Vp ~ — E w ( k ) [aoao4 a -k + aS4 a k«-k] ■ (47b) 

k^O 

Now, if we take the origin of energies at the Gross- 
Pitaevskii value v(0)N(N - 1)/2V (which, it should be 
noted, is an exact eigenvalue of the operator Vh in the 
state ^(iV)) of Eq. (|39| ). then the Hamiltonian can be 
written as a sum of independent contributions from dif- 
ferent values of k: 

^E^ ( 48 ) 

k^O 

where: 

-1/4 + \ V S ) / t t 

#k = 2 £ k( a k a k + a_ k a -kj + -yy- (,aQa a^ak 

+ aQa a^L k a-k + 4 a -k a o a o + 4 a o a k a -k) ■ (49) 

In the next Section, we shall restrict our attention to 
the single- mode Hamiltonian H^. We shall see that the 
ground state of can be found quite easily using a 
number-conserving variational ansatz^ leading to re- 
sults for the ground state energy and excitation spectrum 
that are very similar to those of the standard (number 
non-conserving) Bogoliubov theory. The implications of 
the diagonalization of on the diagonalization of the 
full many-body Boson Hamiltonian H will be discussed 
in more detail in Sec. IVII1 
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V. NUMBER-CONSERVING APPROACH: 
VARIATIONAL FORMULATION OF 
BOGOLIUBOV'S THEORY 

In view of the conceptual difficulties of the number 
non-conserving formulation of Bogoliubov's theory re- 
viewed in Section IIII1 we now want to study an alter- 
native formulation of this theory in which the number 
of bosons is conserved, and which therefore should be 
free of many of the problems encountered in the num- 
ber non-conserving version. This formulation, which has 
already been presented in slightly different form in the 
literature^^ consists in using a variational ansatz to 
find the ground state wavefunction of each single-mode 
Hamiltonian independently from the Hamiltonians 
.ffk'(^k) °f the other momentum modes and writing the 
ground state wavefunction of the total Hamiltonian H as 
a tensor product of all these separate single-mode ground 
states. However, since the previous expositions of this 
method are not widely in use, we here shall present a 
comprehensive review, which will allow us to assess the 
validity of several aspects of the standard, number non- 
conserving Bogoliubov formulation. 



A. Ground state energy of the single-mode 
Hamiltonian 

We now proceed to diagonalize the Hamiltonian 
independently of the remaining contributions ifw^k) to 
the total Hamiltonian H of the system; in other words, 
we consider a fictitious system where bosons are only 
allowed to be in one of the three single particle states 
with wavevector 0, -l-k or — k. We shall therefore restrict 
our attention to trial states of the form (without loss of 
generality, throughout this paper it will be assumed that 
the total number of bosons N is even): 



N/2 



(50) 



71=0 



where each basis state \N — 2n, n, n) has (N — 2n) bosons 
in the condensate, and n bosons in each one of the two 
momentum states ±k. To clarify what we mean exactly 
by the notation \N— 2n, n, n) in the present context, if |0) 
designates the vacuum state for bosons, the normalized 
general state | N — n — m, n, to) can be defined as: 

\N-n-m,nM= i °> ( , k) , |0)- 

(51) 

The expectation value of Hk in the un-normalized state 
\tjfk) of Eq. ([50)) is given by: 



(Huh 



(V'klV'k) 



(52) 



Let us calculate the above expectation value. To simplify 
the notation, in the following we shall denote by |n) the 
ket \N — 2n,n,n). We then can write: 



(53a) 



a lk a k|"> = a'_k a -kl n ) = n \ n )i 
a o a o a ik ak l n ) = a o a o a -k a - k l n ) = — 2n)|n), 

(53b) 

alal_ k a a \n) = {n + l)y/(N - 2n)(N - 2n - Tj\n + 1), 

(53c) 

aka-k4 a ol n ) = n V( N - 2n + 1)(N - 2n + 2)\n - 1), 

(53d) 

Using the above equations, one can easily show, after a 
few manipulations, that: 



v(k) 
~2~V~ 

v(k) 
~2~V~ 



N/2 

C « [ n£ k + ^n(N - - 2// i ;/) 

n=0 



v(k) 



N/2 



2m + l)|m) 



rn—l 
(N/2)-l 



( l + l)Ci+W(N -21- 1)(N - 2010 

(54a) 



Now, it is not difficult to see that the second sum on the 
rhs of the above equation can be extended to to = 0, 
since the factor mC m -\ in the summand will make the 
extra term vanish. Also, the third sum can be extended 
to I = N/2, because the factor (N — 21) inside the square 
root will make this term vanish. This allows us to write: 



ffkl^k) 

v(k) 
~~2~V~'' 

v(k) 

~2~V 



N/2 



n=0 



v(k) 



n(N - 2n) 



+ 



+ 



',C n _i y/(N -2n + 2)(N - 2n + 1) 



(n + l)C n+iy /(N -2n- 1)(N - 2n)}|n). (55) 



Using this last result, we obtain that the expectation 
value (ipk\H k \ipk) is given by: 



v(k) 



n(N - 2ri) 



N/2 

(^iHki^k) = E{i c '"i 2 

n=0 

+ ^-nC*C n _ iy /(N -2n+2)(N-2n+ 1) 

+ ^(n + l)C*C n+iy /(N -2n-l)(N- 2n)}. (56) 

Let us assume, for simplicity, that the coefficients C n are 
real. Then, for v(k) > 0, we see that the expectation 
value ( , 0k|#k|V'k) will be lowered if the coefficients C n 
have alternating positive and negative signs. In this case, 
the terms on the second and third line will be negative, 
making the expectation value lower than what one would 



9 



obtain if consecutive coefficients have the same sign, in 
which case products of the form C n C„_i and C n C n+ i 
will be positive. Below, we will show that Bogoliubov's 
theory corresponds to the following geometric ansatai^ 
for the coefficients C n : 



C n — C (— c k) r 



(57) 



where the constant Ck is to be determined variationally 
(Co will turn out to be an overall constant which cancels 
out in the normalization of and whose value is hence 
unimportant for the evaluation of expectation values of 
physical observables) . Note that, for the C„'s to have 
alternating positive and negative signs, Ck has to be pos- 
itive. On the other hand, we expect on physical grounds 
that the coefficients C n will decrease with increasing val- 
ues of n, or, in other words, that the probability ampli- 
tude of states |n) with a large number n > 1 of bosons 
outside the condensate will be small. This implies that 
the constant Ck must be less than unity. It then follows 
that Ck is subject to the following constraint: 



< c k < 1. 

Inserting the variational ansatz (|57[) into Eq. 
making use of the approximation: 



^N(N+l) = N^l + ±^N+± 



(58) 
56|) . and 

(59) 



which is valid for N 3> 1, we can rewrite Eq. (|56[) in the 
form: 



(V-kl^kl^k) 



N/2 

2 E 

n=0 



v(k) 
~W 
v(k) 

~2V 



v{k) 



{(ck) 2 "[«e k +— »! V .2,/) 
"(ck) 2 "" 1 ^ 



N — 2n 



) 



(n + I)(c k ) 2 



(60) 

The summations in the above equation can be calculated 
analytically by taking successive derivatives with respect 
to the variable x of the finite geometric sum: 



N/2 

H xn = — 

n=0 



hence obtaining the following results: 
f 



N/2 

£'< 

n=l 
N/2 



X 2 



+ 1 



E 



n x = x 



(1-x) 2 

1 — X 2 +± 



N 
~2 



1 



.r 



N/2 



X. 



(1- 

2(1 



N 
2 



1 

x N/2 



(61) 



(62a) 



1 



x 



r l+N/2 



(1 



V 2 



r N/2 



(1 



N / N 



(62b) 



Now, it would be highly impractical to use the full ex- 
pressions on the rhs of Eqs. (|62p in the evaluation of the 
sums on the rhs of Eq. (|60|) . Luckily, these expressions 
simplify considerably for N — > oo and < x < 1, in 
which case we can write: 



N/2 

n=l 
N/2 



(1-x) 2 ' 

El n ^ 
nX ~ (1 - x) 2 ' (1-x) 3 ' 



2x 2 



(63a) 



(63b) 



Using these last two equations in Eq. (|60|) . we obtain 
(we remind the reader that ns = N/V is the density of 
bosons in the system): 



(V'klffklVk) = |C | 2 e k + v(k)n B (l-—)\J2 



N/2 



k ) 



(64a) 



\ C o\ 2 {j r ^[^ + v(k)n B ]-v(k)n Bj] -^} 



(i-4) 2 



(64b) 



On the other hand, from Eq. , we easily see that the 
norm of the wavefunction (V'klV'k) is given by: 



N/2 N/2 



n=0 

Co| 2 



n=0 



i-4' 



(65a) 
(65b) 



where, in going from the first to the second line, use has 
been made of the asymptotic form of Eq. (|61[) . namely: 



N/2 

E^ 



i 



l-x' 



N ->• oo, < x < 1. 



(66) 



Now, if we divide Eq. (|64b|byEq. (f65b|) . we can write 
for the expectation value (-ffk)k of Eq. (|52[) the following 
expression (notice that | C 1 2 cancels out in the division) : 



k/k 



i-4 



fe k + v(k)n B ] - v{k)n B 



Ck 



1 - c 2 ' 

1 L k 



(67) 



Minimization of the above expectation value with respect 
to Ck leads to the following quadratic equation)!^ 



4-2 



v{k)n B 



c k + 1 = 0, 



where we defined 



£k = £k + v(k)n B . 



(68) 



(69) 



The above quadratic equation has the following two 
roots: 



v{k)n B 



± 



fk 



v{k)n B 



- 1. 



(70) 
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FIG. 1: Plot of the constant c k of Eq. (|71[) vs. wavevec- 
tor k. Here, k is the dimensionless combination k = 



hk/ ^/2mv(k)riB ■ 



particle state of momentum k is given by: 

(V'k|a k ak|V'k) 



N k = (a k ak)k 
Now, it is easy to see that: 



(V'klV'k 



(74) 



N/2 N/2 

(V>k|a k a k |Vk> = |Co| 2 £ n \°"\ 2 = \Co\ 2 £ ™4", 

n=0 n=0 

(75a) 

„2 



|Co| 



(75b) 



where, in going from the first to the second line, use has 
been made of Eq. (|63al) . Now, using expression (|65bD . 
we finally obtain: 



Of these two roots, only the one with the minus sign 
obeys the constraint < c k < 1 of Eq. (|58|) for arbitrary 
values of £ k . We shall therefore writes^" 2 ^ 



1 



Ck 



v{k)n B 



3 



v(k) 2 n B 



(71) 



This result for the constant c k , which is plotted as a 
function of the wavevector k in Fig. [TJ fully determines 
the coeffcients C n — Co(— c k ) n of the variational ground 
state wavefunction |^ k ) of the Hamiltonian i? k (apart, of 
course, from the overall constant Co). This, in turn, will 
allow us to determine the expectation values of physical 
observables in the ground state \ipk). In particular, the 
expectation value of iJ k in the ground state |?/> k ) can 
readily be found if we use the result (|71[) for c k in Eq. 
(|67p . upon which we obtain: 

(ffk)k = -^(ek + n B i;(fc)-£ k ) < 0, (72) 



with the definition: 



£k[£k + 2n B v(k)]. 



(73) 



The result ([72]) is exactly the result one obtains in the 
standard, number non-conserving Bogoliubov approach 
for the expectation value of a given contribution _ff k to 
the Bogoliubov ground state energy. We shall discuss 
the meaning of this observation in more detail below. For 
the moment, we want to use our result for the variational 
wavefunction |^ k ) to explore the properties of the ground 
state. 



Variational result for the depletion of the 
condensate 



We now turn our attention to the depletion of the con- 
densate. The average number of bosons -/V k in the single- 



\a k a k ) k = 



i-4' 



(76) 



If we now use the expression of c k in Eq. (J7TJ), we find 
after a few manipulations: 



£k + n B v{k) 



1 



(77) 



This is again the same result one obtains in Bogoliubov's 
standard, number non-conserving approach. Unfortu- 
nately, and as we have already discussed in Sec. IHI Al 
the above expression of Af k diverges for k — ¥ 0, which 
does not make much sense for a system of A~ bosons. It 
is indeed quite easy to see that, A^ k being the ratio of two 
finite sums: 



(78) 



with < Cfe < 1 , the ratio must remain finite for all val- 
ues of Ck in the range (0,1). Luckily, the origin of the 
divergence in Eq. ([77)1 can be easily elucidated, and has 
to do with our use of the approximate limiting expres- 
sions of the sums (|63a|) and (fB"6"j) (which both diverge as 
c k — > 1) when evaluating the average (a k a k ) k . If we use 
the full expressions (|62al) and (|6"T|) instead, we find the 
result: 



A^ k = 



1 



<T + 1 



n 2 + N 



1 



(79) 



It can be verified that the above result, which to the 
best of the author's knowledge has not been studied pre- 
viously, is finite for all values of c k . In particular, as 
k — > 0, c k — > 1, and a Taylor expansion in x — c k near 
x = 1 shows that A^ — > N/4. This behaviour is shown 
in Fig. [5] where we plot the depletion Nk from Eq. (f79")) 
as a function of the wavevector k — |k| for a number of 
values of the total number of bosons N. 



11 



100 
90 
80 
70 
60 
50 
40 
30 
20 
10 








0.02 0.03 

k 



FIG. 2: Depletion Nk = (ffl k £ik} as obtained from the varia- 
tional approach, Eq. (|79[) , for various values of N. The var- 
ious curves are for N = 100, 200, 300 and 400 bosons, from 
bottom to top. As it can be seen, the depletion iV k — > N/4 
as k — > . Here aga in, k is the dimensionless combination 
k = hk/ y/2mv(k)riB- 



The result of Eq. (jT§)) allows us to obtain more mean- 
ingful results for expectation values of one-body opera- 
tors of the form O — O k a\ji k than were obtained 
within the standard, number non-conserving Bogoliubov 
approximation. To revisit the example given in Sec. 
IIII Al let us consider the Fock part of the Hamilto- 

nian Hp ~ \ Sk^o n Q v O s -)( a l. a ^ + a Lk a -k), where for 
simplicity we made use of Bogoliubov's approximation 
ao = aj ss y/No- While in the standard Bogoliubov ap- 
proximation the expectation value of a given mode k of 
Hp in the Bogoliubov ground state diverges like 1/k as 
k — > (which, as we discussed earlier, is problematic), in 
the variational Bogoliubov method this expectation value 
is finite for all values of the wavevector k, and is given by 
the usual expression involving the product of the number 
Nq of bosons in the mode k = by the number iVk of de- 
pleted bosons in each of the modes k and — k (since these 
are the only modes kept in Hp), times the interaction 
energy v(k)/V. The fact that all quantities involved are 
well-defined and finite within the variational approach 
is reassuring, and constitutes a significant improvement 
with respect to the standard, number non-conserving Bo- 
goliubov approximation. 



Variational result for the anomalous averages 

(a k a_ k ) and (a k al k ) 



the specific value of N (in other words, the constant Ck 
is independent of N for A" — > oo). Under these circum- 
stances, we shall write for the wavefunction of a system 
of (N — 2) bosons the following variational ansatz: 

(JV-2)/2 

\MN-2))= ^ C n \(N-2)-2n,n,n). (80) 

71=0 

It is now easy to verify that: 

N/2-1 
n=0 



x \N- 2(n + l),n + l,n + l), 



and hence that: 



N/2 

nC n -i\n), 

n=l 



N/2 



(81a) 
(81b) 



(MN)\alal k \MN - 2)) = £ nC n _iC*, (82a) 



C, 



Ck 



^ (82b) 



Using the limiting form (|63a[) for the geometric sum as 
A" — > oo, we obtain: 



(i> k (N)\aial k \i; k (N-2))^-C 2 - 



(83) 



On the other hand, the limiting form (f6"T|) as A" — > oo of 
the norm of the wavefunctions ^^(N)) and \tpk.(N — 2)) 
is given by: 

(^k(A0lV>k(A0) = (MN - 2)\MN - 2)) = Cl 



(84) 

If we use the normalized ground state wavefunctions 
|^ k (A0) such that: 

|V>k(A0) 



(85) 



then the anomalous average (o^a^)^ can be defined in 
the following way: 

(«k«-k)k = (MN)\atal k \MN - 2)). (86) 
Using Eqs. ((83]) and JM]), we find: 



We now want to investigate the anomalous average 
(aka_k:) m the variational formulation of Bogoliubov's 
method. In order to do so, we note that for A^ ^> 1 the 
coefficients C n of the variational wavefunction IV'k(AO) = 
Y^n=o^n\ n ) do not explicitly depend on the number of 
bosons N, and are given to a very good approximation 
by the ansatz expression C n — Cq(— Ck)™, regardless of 



( a k a -k)k 



Ck 



i-4' 



(87) 



Given the expression (|71[) of Ck, we obtain, after a few 
manipulations: 



( a k a -k)k 



n B f (k) 
2E k ' 



(88) 
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FIG. 3: Anomalous average (a k ffl-k) as obtained from the 
variational approach, Eq. (|91[) . The Bogoliubov result is 
shown as the dotted line, and diverges negatively in the k — > 
limit. The solid line is the result of the variational approach 
for N = 100 bosons, and the dashed line is the variational 
result for N = 200 bosons. As it can be seen, in the variational 
method the quantity (dkO-k) goes to a finite limit as k — > 0, 
as it should. (As in Figs. [T] and [21 k in this figure is the 
dimensionless combination k = hk/ y / 2mv(k)nB ■) 



an expression which again has not been studied previ- 
ously, and which can be shown to be bounded for all 
values of Ck in the range (0, 1) (see Fig. |3]). 

The result of Eq. (J9TJ) allows us again to obtain 
more meaningful results for expectation values of one- 
body operators of the form O = Ok^ka-k or 
O = Oka k a! k than were obtained within the stan- 
dard, number non-conserving Bogoliubov approxima- 
tion. To revisit the example given in Sec. IIII A[ let 
us consider the pairing part of the Hamiltonian Hp 
\ X^k^o n of(k)(a k a^ k + a_k£ik), where for simplicity we 
again made use of Bogoliubov's approximation ao = ~ 
VNq. While in the standard Bogoliubov approximation 
the expectation value of a given mode k of Hp in the 

Bogoliubov ground state is given by -[n B v(k)Y/(2E k ) 
and diverges like — 1/k as k — > 0, in the variational Bo- 
goliubov method this expectation value is finite for all 
values of the wavevector k, and this even in lower spatial 
dimensions. 



D. Elementary excitations: variational approach 
for the single-mode Hamiltonian i/ k 



which is nothing but the standard Bogoliubov result of 
Eq. (|23|) above, which diverges negatively like — 1/k as 
fc^O. 

We now show that the divergence of the rhs of Eq. ([55)1 
as k — ► is also a direct consequence of using asymptotic 
forms as N — ¥ oo of geometric sums, and that, if the 
finite N results are used throughout, then a completely 
diffrcnt result will be obtained which remains finite for 
all values of the wavevector k. Indeed, using the notation 
x = c k , we can write: 

(Vk(A0|a k aL k |Vk(iV-2)) 



r 2 

Ck 



E 



nx 



r 2 

Ck ' 



1 - x 1+N ' 2 



{1-xf 



N 
~2 



n=l 
T N/2 



1 



(89) 



In a similar fashion, we obtain for the norm of the wave- 
functions \4>k(N)) and \ipk(N — 2)) the following expres- 
sions: 



<Vk(A0|Vk(A0) = cl 
(MN-2)\MN-2)) = Cf ) 



1 _ X 1 + N / 2 
1 -x '■ 
1 - x N / 2 
1-x 



(90a) 
(90b) 



Combining Eqs. (j8TJ)) and (|9T)]) allows us to write the 
final expression: 



Ck 



(l-c^Xl-c^ 2 ) 



1 - c 



N+2 



N 



(91) 



We now introduce the following operator : ?? 



Q k = tt k akao + «ka aL k - 



(92) 



We want to choose the constants tik and "Dk in such a way 
that acting on with a k creates an excited state of 
the Hamiltonian H^. A necessary condition for that to 
happen is that a£|'0k) has to be orthogonal to |?/>k), i-e., 
that: 



(^ k |(4l^k)) 



= 0. 



(93) 



Taking the adjoint of the above equation, we obtain the 
condition 

(tfklaklV'k) = 0, (94) 

which will always be satisfied if we require that 

akIV'k) = 0. (95) 

Let us calculate the effect of the action of au on |^k)- 
We have: 

N/2 



«k|^k> = Y. Cn [v-ky/n(N- 2n+ 1) 

n=0 

x \N - 2n+ l,n- l,n) 



+ v k y/(n + l)(N -2n)\N - 2n - l,n,n+ 1) 

(96) 

Noting that the n = term in the first summation and 
the n = N/2 term in the second summation vanish, we 
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relabel n n — 1 in the second sum, hence obtaining: 

N/2 



aklVk} = £ [c„iikVn(iV-2n+l) 



+ C„_iS kV / n(iV - 2n + 2) |JV - 2n + 1, n - 1, n). 

(97) 

Requiring that aiklV'k) = leads to the following condi- 
tion: 



Vk N-2n + 2 



u k V N - 2n + 1 



(98) 



For AT 3> 1, the square root is very close to unity for 
most values of n such that 1 < n < N/2, and so we 
obtain that, to a very good approximation, the condition 
CKklV'k) = is equivalent to: 



Vk 



(99) 



Using the fact that C n = Cq{— c k ) n , we finally obtain: 

(100) 



— = c k . 



Hence, we obtain that the operator annihilates the 
variational ground state \ipk) of the Hamiltonian Hk if 
the constants Vk and Uk have the ratio specified in Eq. 
(|100|) . This result is quite important, given that the equa- 
tion a k |-0k) = in effect ensures that the state ajjV'k) is 
orthogonal to IV'k), and is hence an excited state as we 
discussed earlier in this Subsection. 

We now want to establish that the operator a k creates 
an excited state of momentum k. To this end, let us 
find the average momentum in the state ajjV'k), which 
is defined as: 



(p)r 



(V'klakPaklV'k) 



(101) 



with P = ftka k ak the total momentum operator. 
Noting that the commutator [P,a k ] is given by: 

[P,4] =7*4> ( 102 ) 
and that P|"0k) = 0, we can easily write: 

(^klakPttkl^k) = ftk^klakttklV'k), 
from which we finally obtain: 

= Kk, 



(103) 



(p>r 



(104) 



thus proving that the state ajj^k) has an average mo- 
mentum hk. 

Let us recapitulate what we have done so far. We have 
defined a new operator ctk such that the state ajjV'k) is 



(i) orthogonal to \ipk) and (ii) has an average momen- 
tum of hk. These two properties allow us to identify the 
state a^l^k) as an an excited state of the Hamiltonian 
Hk with momentum hk. At this point, we make the ob- 
servation that, while the requirement ak|"0k) = allowed 
us to determine the ratio of Vk to iik, Eq. (jlOOp . these 
constants are otherwise still arbitrary. In order to deter- 
mine their actual values, we shall require that the excited 
state a k |^k) be normalized to unity, i.e. that: 



(^klakOklV'k) = 1, 



(105) 



where we remind the reader that \ipk) denotes the nor- 
malized ground state of the single-mode Hamiltonian Hk- 
Using the fact that = [c*k> 4] + a k a kj the condi- 

tion (|105|) becomes: 



(V'k|[ak,a k ]|V'k) = 1, 



(106) 



where we used the fact that aklV'k) = 0. An immediate 
way to satisfy this last condition is to require that the 
commutator [ak,ct k ] be equal to unity: 



[«k, 4] 



1. 



(107) 



We thus see that the canonical commutation relation 
[ak,a k ] = 1, which is imposed a priori in the standard 
textbook formulation of Bogoliubov's theory, emerges as 
a natural requirement that we need to satisfy in order for 
the excited state ajji/'k) to be normalized to unity. 

Now, if we calculate the commutator [ak, a k ], given 
the definition (f92|) . we find, after a few manipulations: 



[a k , a k ] = (tt k - Wk) a o a o - "k a ik a k + Wk«-k a -k- (108) 

While in Bogoliubov's theory the commutator [ak, cn k ] is 
a c-number, the fact that the rhs of the above equation is 
an operator makes it impossible to satisfy the constraint 
[c*k)C* k ] = 1 exactly in the number- conserving formal- 
ism. However, we can try to satisfy this constraint in an 
averaged sense in the ground state \ipk) by imposing the 
condition given in Eq. (|106|) . which can be rewritten in 
the form: 



([ak, a k ]) k = (u k - *>k) ( a o a o)k - (a k a k)k 



(109) 



where we denote by (•••) k the quantum average 
("0k | ■ • • \ipk)} and where we used the fact that (a k <2k)k = 
( a -k a -k)k- Now, the normalization condition Q106P 
leads to the following constraint on the constants Uk and 

Vk- 

K-vl = 7k, (HO) 

where, for compactness of notation, we defined the quan- 
tity 7k such that: 

1 



ll = 



{a\a )k - (a k °k)k 



t \ ' 
aoa )k 



(111a) 
(111b) 
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In going from the first to the second line of Eq. (Illip we 
assumed that (a k o,k)k "C (ajao)ki that the depletion 
of the condensate into the single-particle state with mo- 
mentum k is small. Note that this is not always the case, 
as in the standard formulation of Bogoliubov's theory the 
quantity Nk = (a k <2k)k diverges in the k — ¥ limit. This, 
incidentally, gives us another reason why it is important 
to be able to formulate a version of the theory in which 
the expectation value (a^ak) always remains finite. As 
mentioned above, in the number non-conserving formu- 
lation of Bogoliubov's theory, (a^ak) diverges as k —> 0, 
and taking this result at face value would lead to a com- 
pletely different result for 7 k in Eq. (jlllbj) . which would 
in turn affect the evaluation of the excitation energies 
given below, Eq. (|116[) . 

Using Eqs. (|100[) and (|110[) , we now can determine the 
values of the constants Uk and Vk- These can be written 
in the form: 



Uk = 7kUk, v k = 7k«k, 

1 c k 



Using the result (|71|) for Ck, we finally obtain: 

2 1 /e k + n B v(k) 
Uj, = - h 1 



2 V 



1 /£ k + n B v(k) 

2 V 



£k 



- 1 



(112a) 
(112b) 

(113a) 
(113b) 



Note that these are the same expressions for the quanti- 
ties itk and Vk defined in the conventional formulation of 
Bogoliubov's theory. The above results, Eqs. (|1 12[) and 
(|113|) . allow us to rewrite the expression of ak, Eq. (|92|) . 
in the form: 



ak = 7k I w k a k aj + v k aoa'_ k 



(114) 



If in the above expression we replace do and a by y/No, 
and then use the fact that 7k — 1/ V^o as indicated by 
Eq. (jlllbj) . we recover the usual expression, Eq. (O, of 
ctk in terms of ctk and al k used in the standard, number 
non-conserving formulation of Bogoliubov's method. 

We now are finally in a position to find the energy of 
the excited state a k , which is given by the expression: 

AT?Wn\ WgkgkQkljk) /7|rrl7\ fit* \ 

& E excW = —rT, T7~~ x — (W-HklW, ( 115a ) 

(Vklakakl^k) 

= $k|ak#k4l&) - (^kl^kl^k)- (H5b) 
Here, in going from the first to the second line, we used 
the fact that the excited state a k l^k) is normalized to 
unity, i.e. that (^klak^kl^k) = 1- In Appendix lAl we 
show that, for weak perturbation, the Hamiltonian Hk 
can be written in the form: 29 

H k ~ ~2 [ £ k + n Bv(k) - E k ] + ^E k (ala k + aL k a -k), 



where Ek is the Bogoliubov spectrum given in Eq. (|73l) . 
This last result, when used in conjunction with Eq. (|115p 
and the commutation relations for the ctk's, Eqs. (|A17jl - 
(|A19[) of Appendix [S] leads to the conclusion that the 
excitation energy of the Hamiltonian Hk from the ground 
state \4>\s) to the excited state a k |7/> k ) of momentum +k 
is given by the following result: 

A£« (k) = ~£ k = ye k [e k + 2n B v(k)]. (117) 

We here would like to emphasize that this is the excita- 
tion energy of the Hamiltonian Hk . The full Hamiltonian 
H = X)k^o ^ k contains identical contributions from Hk 
and -ff-k 7 and hence one can easily see that the total 
energy cost A2?*£*(k) of the excitation a^\ipk), account- 
ing for contributions from Hk as well as -ff_k, is twice 
the amount quoted in Eq. (|117[) . and is given by the 
standard Bogoliubov result, namely: 

AE% c (k) = A£<2(k) + AEgl (-k) = E k . (118) 

Before we end this Section, we briefly comment on why 
we define the excitation operator ak by the expression 
(|92p . and why we find it useful to impose the condition 
(|95|) as a way to fix the value of the constants Uk and 
Vk. As a matter of fact, it is not difficult to verify that 
the operator a^ao also creates a state which is orthogonal 
to the ground state \ifik) of the operator Hk, and could 
therefore be a possible choice for an excited state of the 
Hamiltonian Hk- Indeed, the action of the operator a k ao 
on gives: 

N/2 

4«o|^k) = C Wi N - 2n )(n + l)|JV-2n-l, n+1, n). 

n=0 

(119) 

Now, the rhs of the above equation belongs to the Hilbert 
space spanned by kets of the form \N — 2n— 1, ra + 1, ra), 
which is totally disjoint from the Hilbert space spanned 
by states of the form \N — In, n, n) and to which the 
ground state belongs, hence the orthogonality of the state 
a k ao|'0k) and the ground state ji/'k) of the Hamiltonian 
Hk- By the same token, (V'k|cfc k J'0k is automatically zero, 
and so ajji/'k) is automatically orthogonal to \if?k). So, 
why do we choose the special combination ctk = Uk&kflo + 
v k£t_k a o instead of simply ak = &k&o? And why do we 

take the trouble of requesting that ak\ipk) itself has to 
vanish? 

To answer the first question, we note that the oper- 
ator a k = Mka k a o + "Ska-k^o increases the momentum 
of the system by an amount +fik. The expression of 
a k simply emphasizes the fact that this can be done in 
two different ways. The first way consists in removing 
a boson from the condensate and adding a boson in the 
single-particle state of wavevector k, hence the UkO^ao 
(116) contribution to a k . The second way consists in removing 
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a boson of wavevector — k from the system and adding 
a boson to the condensate, hence the term v^a^a in 
the expression of a*. The use of the special combination 

Uk a k a o + ^kG-kGo is convenient in that it allows us to 
combine the two types of processes into one excitation op- 
erator (it also allows us to choose values for the constants 
Uk and -Dk such that the ket ajj^k) is automatically nor- 
malized to unity) . As for the second question, the reason 
we impose the condition (|95[) is in fact to ensure that all 
states obtained by successive application of the operator 
a k are orthogonal to the ground state. Indeed, without 
the condition (|95|) . there is no guarantee that a state of 
the form (a^) 2 !^) satisfies the orthogonality condition 
(i ; k\( a L) 2 \''pk) = 0, even though the state ajjVk) satisfies 
( - 0k|akl , /'k) = (V'klcKklV'k) = 0. Equation (J95J) is useful 
in that it ensures that all states of the form (al) n \tp'k) 
are orthogonal to the ground state iV'k) of the Hamil- 
tonian H^, a fact that can be easily verified by taking 

the adjoint of the expression (V'k|( Q; k)"l'^' k )- fa can m 
fact be shown that the condition (1551) ensures that states 



of the form (a k ) n |V>k) and (a£) m |V'k) are orthogonal to 
each other if n =/= m. The proof, which is quite straight- 
forward, will not be presented here.) 



VI. EXACT NUMERICAL DIAGONALIZATION 
OF THE SINGLE-MODE HAMILTONIAN H k 

Having discussed a number-conserving variational ap- 
proach to the single- mode Hamiltonian of Eq. (j49| . 
we now proceed to perform an exact numerical diagonal- 
ization of this Hamiltonian. Such a numerical treatment, 
which to the best of the author's knowledge has not been 
attempted before, will allow us to assess the validity of 
the variational treatment of the previous Section, and 
possibly point to areas where the variational approach 
may need improvement. We shall again use the com- 
plete basis \n) = \N — 2n,n,n) formed by states with 
(N — 2n) bosons in the k = state, and n bosons in each 
of the states k and — k. Since n can only take values 
between and N/2, the basis {\n}} has (1 + N/2) kets, 
and hence diagonalizing requires the diagonalization 
of an (1 + N/2) x (1 + N/2) matrix. For a system of 1000 
bosons, this would be a 501 x 501 matrix, which is not 
very large given nowadays computing capabilities. 

Using Eqs. (|53p . we can easily write the following ma- 
trix elements: 



m\Hk\n 
v(k) 



vfa) 
V 



n(N - 2n) 



+ 



2V 



(n + l)y/(N-2n)(N-2n-l)S min+1 



+ ny/{N -2n + 1)(N -2n + 2)8 m ,n-\ 



(120) 



We now introduce dimensionless units, where we mea- 
sure energies in units of v(0)ub = v(0)N/V. Hence, we 



shall write: 
(m\Hk\n) 



v(0)ub 

+ lm 



nii< + v(k)—(N-2n) 



(n + 1) 



N 



y/(N-2n){N-2n-l)6 m , n+1 

(121) 



+ jfV( N -2n + l)(N-2n + 2)<5 m ,„_ 

where we denote by £k and u(k) the dimensionless quan- 
tities: 



£k 



S(k) 



?k 



v(0)tib ' 

_ vM 

v(0)- 



(122a) 
(122b) 



Having in mind an interaction potential of the form 
v(r) = gS(r), whereby u(k) = g, throughout this Sec- 
tion we shall present numerical results using for u(k) the 
value 



v(k) 



i. 



(123) 



Now, using the dimensionless matrix elements in Eq. 
(|12ip . it is a relatively easy task to diagonalize the Hamil- 
tonian //k for a given value of the total number of bosons 
N and for various values of the dimensionless variable 
ek- In the following, we shall present plots of physical 
quantities of interest as a function of the dimensionless 
wavevector k such that: 



where we defined: 



k 



fc 



y / 2mnBv(0) 



(124) 



(125) 



A. Ground state energy 

The first quantity we shall be interested in is the 
ground state energy of the Hamiltonian H^. Numeri- 
cally, the ground state energy Eq(N) of the single-mode 
Hamiltonian is defined as the smallest eigenvalue of 
i?k- In Bogoliubov's method, the ground state energy of 
-ffk is the expectation value of this last quantity in the 
normalized ground state |"0k), and is given by Eq. ((721) . 
Fig. 2] shows plots of Eq(N) obtained by numerical di- 
agonalization of i?k for N — 200 bosons on one hand, 
and by using the Bogoliubov expression for (flk)k, Eq. 
(fT2")) . on the other. The agreement between Bogoliubov's 
method and the exact numerical treatment is quite im- 
pressive. The Bogoliubov results are practically indistin- 
guishable from the exact numerical ones for all values of 
fc, except in a narrow interval near fc = 0, where there 
is a small deviation between the two results. Here, we 
shall not spend too much time trying to understand this 
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FIG. 4: Upper panel: Ground state energy of the Hamilto- 
nian for TV = 200 bosons as a function of the wavevector 
k. The circles represent the exact result obtained by numer- 
ical diagonalization of the Hamiltonian H^, while the solid 
line is the result of the standard and of the variational Bo- 
goliubov formulations. Lower panel: Detail from the above 
figure at small wavevectors showing that the exact numerical 
result (solid line) very slightly deviates from the results of 
Bogoliubov's theory (dashed line) near k — 0. 



(rather small) deviation. Suffice it to say that it is near 
k = where Ck — » 1 that we previously found that the 
standard Bogoliubov theory breaks down, leading to un- 
physically diverging results for the averages (a^a^k and 
(akfl-k)k- Here, we speculate that the deviation of Bo- 
goliubov's result from the exact one for E%(N) may be 
due to similar reasons; in other words, this small devia- 
tion may be due to the ground state energy being evalu- 
ated using the truncated version of the sums (|63a[) - (l63bj) 
and (po) instead of the complete ones in Eqs. (|62al) - (|62b|) 
and (EB. 



Ground state wavefunction and depletion of the 
condensate 



While Bogoliubov's theory gives results for the ground 
state energy Eq(N) that are in excellent agreement with 
the exact diagonalization of the single-mode Hamiltonian 
i?k, there is only limited agreement between the exact re- 



FIG. 5: Upper panel: Coefficients C n of the normalized 
ground state wavefunction \ipk) of the Hamiltonian vs. 
the index n for TV = 20 bosons and k = 0.1. The circles are 
the results of the exact numerical diagonalization, and the 
stars are the result of the variational approach. The triangles 
show the ratio C n +i/C n vs. n. Contrarily to the variational 
formulation of Bogoliubov's theory where this ratio is found 
to be a constant given by Ck (see the constant line of squares 
in the figure), in the exact treatment the ratio C n +i/C n is an 
increasing function of n and goes to zero asn-> N/2. Lower 
panel: Same as the upper panel for N = 200 bosons. 



suits for the properties of the ground state wavefunction 
itself and the results of Bogoliubov's method. This is 
a well-known feature of variational methods in general, 
which may, with appropriate choice of trial wavcfunc- 
tions, give excellent approximations to the ground state 
energy of the system, but quite often only produce ten- 
tative agreement as far as the spatial properties of the 
wavefunction itself are concerned. In the case at hand, 
and as it can be seen in Fig. O we find that the coeffi- 
cients C n of the normalized ground state function IV'k), 
as obtained by exact diagonalization of H^, do alternate 
in sign between positive and negative values as predicted 
by the variational formulation of Bogoliubov's method. 
However, the ratio C n +i/C n does not assume a con- 
stant value as the latter theory predicts, which indicates 
that Bogoliubov's method produces good qualitative but 
only approximate quantitative agreement with the exact 
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FIG. 6: Comparison between the depletion of the conden- 
sate iVk = (rajj-aOk vs. wavevector k for N = 200 bosons, 
as obtained from the exact numerical diagonalization of 
(circles), the variational number-conserving formulation of 
Bogoliubov's method (stars), and the standard number non- 
conserving Bogoliubov theory (solid line). 

treatment as far as the ground state wavefunction |^k) 
itself is concerned. More specifically, we find that the ra- 
tio Cn+i/Cn goes to zero at large values of the index n in 
the numerical treatment, which indicates that in the ex- 
act solution the highly depleted states \N — 2n, n, n) with 
n ~ N/2 are much more strongly suppressed than in the 
variational Bogoliubov method. This in turn will have a 
pronounced effect on the depletion of the condensate. 

In Fig. |5l we plot the depletion N k vs. wavevector k 
that is obtained for N — 200 bosons using three different 
methods. The circles show the exact result: 

N/2 

Wk = <4«k)k = ^n|C ; n | 2 . (126) 

n=l 

where the C ra 's are obtained by exact numerical diago- 
nalization of H^. The stars on the other hand, show the 
result obtained from the number-conserving variational 
formulation of Bogoliubov's method, Eq. (f79|) . while the 
solid line shows the result of the standard (number non- 
conserving) Bogoliubov treatment, Eq. ([77)l . It is seen 
that, while the latter result diverges near k — > and is 
therefore unphysical, the variational result is finite for all 
values of k, and goes to the limiting value N/4 = 50 near 
k = 0. It is also seen that the variational method over- 
estimates the depletion of the condensate near k = by 
a factor of about 7 for the particular value of N chosen. 
As we mentioned earlier in this paragraph, this is due to 
the overestimation in Bogoliubov's variational treatment 
of the coefficients C n for large n, which also contribute 
to the calculation of the depletion of the condensate ac- 
cording to Eq. (fT2"6) . 

In Fig. [7J we plot the depletion that is obtained by 
exact numerical diagonalization of the Hamiltonian H k 
for a number of values of the total number of bosons 
N. Again, comparison with the results shown in Fig. [5] 



shows that the variational method consistently overesti- 
mates the depletion near k = for all the values of N 
considered by about one order of magnitude. The deple- 
tion obtained by exact numerical diagonalization of the 
Hamiltonian is seen to be extremely small, and is in 
fact much smaller than what the standard Bogoliubov 
theory predicts in the |k| — > limit. 
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FIG. 7: Depletion of the condensate = (a^a^k as ob- 
tained from the exact numerical diagonalization of the Hamil- 
tonian H k for N = 400, 300, 200 and 100 bosons, from top to 
bottom. 



C. Elementary excitations 

The exact numerical diagonalization of described in 
the two previous Subsections gives access not only to the 
ground state energy, but also to all excited eigenvalues 
and eigenfunctions of the system in the Hilbert space 
spanned by the kets \n) = \N — 2n, n, n), where < 
n < N/2. Such excited states can typically be generated 
by successive application of the operator a k al k (to be 
specific, here the a k 's are the operators defined in the 
number-conserving version of Bogoliubov's theory, Eq. 
([52])). Indeed, if we calculate the quantity a£aL k |'0] £ ), 

where is of the form |^k) = J2n=o C n \N ~ 2n, n, n), 
we find successively: 

N/2 

"klV'k) = An \ N ~ 2n + 1, ™, n-l), (127a) 

71=1 

A n = C n -xu^n{N - 2n + 2) + C n v k ^n(N - 2n + 1), 

(127b) 

then: 

N/2 

a -k a kl^k) =j2B n \N-2n,n,n), (128a) 

n=0 

B n = A n Uky/n(N- 2n+ 1) 
+ A n+1 v ll y/(n + l)(N-2n), (128b) 
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where, in the last equation, it is understood that 
Ai+n/2 = 0- Eq. (|128a[) shows that the doubly excited 

state aj c aL k | , 0k) belongs to the Hilbert space spanned 
by the kets |n) = — 2n,n,n), and hence is accessible 
through the numerical diagonalization of the Hamilto- 
nian H^. 

Let us step back for a moment and try to calculate the 
analytical expression of the energy required to excite the 
system from the ground state |-0k) to the doubly excited 
This is the quantity: 



state a k at]JV , k 



A£E(k) 



(V>k|aka-k#ka k aL k |^k) 
(-0 k |aka~ka k aL k |V'k) 



(^klffklV'k) 



(129) 

Using the diagonalized expression of in terms of the 
ctk's, Eq. (|116p . we find, after a straightforward calcula- 
tion (in which the approximation [ckkiQ^] — 7k a o a o ~ L 
see Eq. (|A19I) from Appendix [2] is used throughout) 

that AEexc{k) is given by: 



A£« (k) = £ k . 



(130) 



This is not at all surprising. Indeed, since the excita- 
tion energy for a singly excited state with respect to 
the Hamiltonian was found to be E-^/2, here for 
a doubly excited state we expect an excitation energy 
2 x (Ek/2) =E k . 

We now go back to numerics, and show that the energy 
of the first excited state in our numerical treatment co- 
incides with the energy of the double excitation given in 
Eq. (|130[) . As we mentioned in the opening paragraph of 
this Subsection, dealing with doubly excited states of the 
form a k al k |i/' k ) has the advantage of keeping us inside 
the same Hilbert space we used to diagonalize the Hamil- 
tonian iJ k . Hence, in our numerical diagonalization pro- 
cedure, if we extract the ground state energy e£\n) 
and the energy of the first excited state E^ (N), then the 

difference AE e ( E(k) = e£\n) - E^ ] (N) should corre- 
spond to the energy of a doubly excited state of the form 
a k a -klV'k) and should be given by the rhs of Eq. (1130[) . 
The upper panel of Fig. [8] shows a comparison between 
the Bogoliubov result (|130|) and the excitation energy 

(2) 

A£ea C (k) as a function of k for N = 200 bosons. The 
agreement is again excellent for most values of k, except 
for small k where the lower panel of that same figure re- 
veals a deviation from linear behavior and the existence 
of a small energy gap near k = 0, where the excitation 
energy goes to a finite value as k — > 0. To probe whether 
this gap, which it should be noted is much smaller than 
the gap predicted by Hartree-Fock theory, is due to a 
finite-size effect, in Fig. |H]we plot the excitation energy 
obtained numerically at k = 0.001 for increasing values of 
the total number of bosons iV" in the system. The results 
obtained show that the value of the gap monotonically 
decreases as the number N of bosons increases, which in 
turn suggests that the gap will go to zero as A" grows 
infinitely large in the thermodynamic limit. 
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FIG. 8: Upper panel: Plot of the excitation energy AE^ c (k) 
vs. wavevector k for TV = 200 bosons. The circles are the 
exact diagonalization result, and the solid line is the result of 
the variational Bogoliubov method for the Hamiltonian H^. 
The dashed line indicates the asymptotic form of Bogoliubov's 
result, v2fc in dimensionless units, as k — > 0. Lower panel: 
Detail of the region near k = 0, where AE^{(k) goes to a 
finite limiting value as k — > 0, signaling the existence of a 
small energy gap at low momenta, which is most likely due 
to finite size effects. 



VII. GENERALIZATION: FULL BOGOLIUBOV 
HAMILTONIAN 

Having studied in detail the number-conserving varia- 
tional formulation of Bogoliubov's method for the single- 
mode Hamiltonian i/k, we now are in a position to 
tackle the more involved case of the full Hamiltonian 
H = J2k^o 01 the interacting Bose system. However, 
before we do so, we want to examine the relashionship be- 
tween what we did so far and previously published work 
on this system. This will be done next. 



A. Connection to previous work 

Up til now, we have been discussing the variational ap- 
proach to the Hamiltonian _ff k describing the interaction 
of the condensate with states of momentum k and — k. 
At this point, we will observe that it is quite remarkable 
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FIG. 9: Plot of the excitation energy AE { ex ' c (k) at k = 0.001 
obtained from the exact numerical diagonalization of the 
single-mode Hamiltonian as a function of the number of 
bosons TV. The fact that the excitation energy decreases with 
increasing values of TV suggests that the small gap seen in the 
lower panel of Fig. [8] is due to finite-size effects and vanishes 
in the thermodynamic limit TV — > oo. 



that most of the results we derived by diagonalizing 
perfectly coincide with the results of the standard Bo- 
goliubov theory, implying that the latter is effectively a 
theory in which each one of the Hamiltonians i?k which 
contribute to the total Hamiltonian H = X)k^o ^k * s 
diagonalized independently from the other Hamiltonians 
-^k'(^k) m essentially disjoint Hilbert spaces. For exam- 
ple, the ground state energy of the system in Bogoliubov's 
theory is given by (note that we still take the origin of en- 
ergies to be the Gross-Pitaevskii value N(N— l)v(0)/2V): 



E 



Bog 



4e< 



£k + n B v(k) - 23 k ) 



(131) 



k^O 



and it can be verified that this quantity is nothing but: 
E Bog = '£$ k {N)\H )c \$ k (N)), (132) 

k^O 

where we remind the reader that \ipk(N)) is the normal- 
ized ground state of the single-mode Hamiltonian iJ k . 
This is rather surprizing, since we expect the ground state 
energy of the system to be given by an expression of the 
form: 

E Boa = (*{N)\H\*(N)), (133a) 
= J2(*(N)\H k \y(N)), (133b) 

k^O 

where | \I/ (TV)) is the ground state wavefunction of the 
total Hamiltonian H . In what follows, we want to exam- 
ine how, within the standard formulation of Bogoliubov's 
theory, one can go from Eqs. (|133|) to the result in Eq. 
(I132p , which will entail uncovering how the ground state 
wavefunction ^(TV)) for the total Hamiltonian H has 



been constructed in the literature from the ground state 
functions \ip^(N)) of the single- mode Hamiltonians H^. 

The main argument allowing one to construct the 
ground state |<J/(TV)) of the total Hamiltonian H from 
the single-mode ground state wavefunctions |^ k (TV)) is 
based on the following reasoning. As we have already 
seen in Sec. [Hi in the standard formulation of Bogoli- 
ubov's method the creation and annihilation operators of 
bosons in the condensate, a J and a , are replaced by the 
c-number \/Nq ~ \/iV. After this replacement is made, 
the various Hilbert spaces where the single-mode Hamil- 
tonians -ffk act are totally disjoint from each other. This 
in turn implies that the single-mode Hamiltonians {H^} 
can be diagonalized independently from one another, and 
that the ground state wavefunction of the system can be 
written as the product of the ground state wavefunctions 
of each of the single-mode Hamiltonians H^: 



|*(iv))=nhMA0>. 



(134) 



Now, it can be shown that, even though the operators 
oq and a have been replaced by y/N in the expression 
of the Hamiltonian H , minimization of the expectation 
value of the single-mode Hamiltonian in the single- 
mode wavefunction \ifa(N)) ~ ^/l — c k X)n=o( — c k)™l n ) 
leads to the same variational coefficients Ck as in Eq. 
([7T|) . This leads, in the notation of the classic work of 
Lee, Huang and Yang£ to the following expression of the 
normalized ground state |$(JV)) of the Hamiltonian H 
(see also Ref. [30h : 



J1 1= 7100=0 

x |ni,ni; . . . jfioo,^), 

where: 



(135) 



\m,m; . 



17 1 % K % |0). (136) 



In the above equations, ni is the number of bosons 
in the states ±k^, and the coefficient C nii .... noo = 
(~ Cki)" 1 (— c k 2 )" 2 ' ' ' (~ Ckoo)™ 00 represents the probabil- 
ity amplitude of finding the system in the many-body 
state |7ii,ni; • ■ • ; jt-oo, "oo). In principle, we would like 
the coefficients C ni ,..., naa of the wavefunction |\l/(iV)) to 
vanish if m + ri2 + ■ ■ ■ + ^oo > N/2: 



if m + n 2 + h noo > N/2. (137) 



However, enforcing the above constraint is not expected 
to give rise to any significant change in the value of the 
ground state energy in the thermodynamic N — > oo limit. 
Going back to Eq. (|135[) . the quantity Z defined in that 
equation is the normalization constant £ 



n 



(138) 
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With the definitions (jl35[) and (|138[) , one can verify that 
the ground state energy E (N) = (^(N)\H\^(N)) is 
given by the rhs of Eq. (|131|) . which is again the same 
as the rhs of Eq. (|132[) . One can also verify that the de- 
pletion of the condensate is still given by an expression 
of the form (176p . namely: 



(139) 



<*(j\oi4o k iw) = 



Finally, if we forego the number conserving expression 
of the operators «k and a k derived in Sec. IV Dl above, 
and use the number non-conserving approximations of 
Eq. ©, 



o=k = u k 4 + t'ka-k, (140) 



with Uk and the quantities defined in Eq. (|15|) , one can 
easily verify that the action of the operator a?k on |$(JV)) 
yields zero, and that the excitation spectrum generated 
by the a k operator is indeed given by the Bogoliubov 

spectrum = \J £k(£k + 2nsvQs)). 

At this juncture, we would like to attract the reader's 
attention to the fact that in the expression (|135j) of the 
wavefunction |\l/(iV)), the number of bosons in the k = 
is not explicitly specified (this is indeed how Lee et al. 
write the ground state wavefunction in their work, Ref. 
7). One may correctly observe that the number of bosons 
in the k = state needs not be written explicitly, be- 
cause once we specify how many bosons ni are in each 
state with wavevector ±k i; then the number of con- 
densed bosons is automatically known and is given by 
[N — 2Y\ rij). But in fact the reason Lee et al. do not 
specify the number of bosons in the k = state is a more 
trivial one, and has to do with the fact that these authors 
use Bogoliubov's approximation of replacing the opera- 
tors ao and aj by the c-number y/N (one place where 
this is made obvious is the expressions of the excitation 
operators they use, «k and a t , which does not contain 
any trace of the operators ao and a , unlike to the ctk's 
used by Leggett in Ref. [2!| which are similar to those 
we used in Sec. IV Dl above), and hence their calculation 
is one where the k = state is removed altogether from 
the Hilbert space used to describe the system. A major 
question that therefore arises is to know if and how the 
standard results that were derived using the expression 
(|135p of |* (AO) wm change if we put the k = back into 
the Hilbert space, and we keep a more accurate tally of 
how many bosons are present in this k — state, hence 
enforcing the conservation of boson number. This will be 
the subject of the next Subsection. 



B. Variational treatment of the full Hamiltonian 

H = Ek^o H * 

We now want to generalize the variational approach 
of Sec. [V] to treat the full Hamiltonian H — X)k^o ^ k 



of the interacting Bose system. To this end, we shall 
use for |\ff(iV)) an expression of the form (compare with 
Eq. (|3"9")l of the single-mode theory; we here denote by 
M the total number of momentum modes kept in the 
calculation, which will eventually be sent to infinity): 



MN)) = z E ••• E c ni c n2 ...c nM 

m— tim=0 
M 

x \N - 2^n i ;n 1 ,n 1 ; . . .;n M ,n M ), (141) 



where the normalized basis wavefunctions are given by 
(compare with Eq. (11361) ): 



M 



\N - 2y j n i ;n 1 ,nx\ ■ ■ .;n M ,n M ) 



i=l 



(142) 



£ (4 )" ! ( flt k T 



2 = 1 



Note that the ground state wavefunction in Eq. (|141[) 
is not a simple product of ground state wavefunctions of 
the single- mode Hamiltonians fZk, as in Eq. (|134l) . and 
that, even though the i/k's may seem to be decoupled, 
the presence of all the rij's in the number of condensed 
bosons [N — 2 ^2iLi n i] acts like an implicit and rather 
nontrivial coupling between all these Hamiltonians. Note 
also that the summations over the rij's in Eq. (| 141[) ex- 
tend from to a value rii jmax , emphasizing the constraint 
that the number of bosons in the k = state in each of 
the basis wavefunctions has to be greater or equal to zero, 

l - e - T,™i n i < N / 2 ( see E q- (USB)- A possible (but, by 
no means, unique) choice for the values max is given 
by: 



TlM.ma; 



= N/2, 
= N/2- 



N/2 - m 



n M -i- 



(143a) 
(143b) 

(143c) 



It can be verified that the above parametrization of the 
ni mal 's exhausts all possible states compatible with the 
constraint (|137[) . In what follows, however, we will not 
find it important to keep track of the constraints (|143p . 
and in the evaluation of physical expectation values we 
shall for simplicity extend the summations to infinity, 
as it is done in the standard Bogoliubov theory of Eq. 
(|135[) . However, we shall find it essential to keep track of 
the number of bosons in the k = state, as this signif- 
icantly modifies some important results of the standard 
Bogoliubov approach, changing the nature of the exci- 
tation sprectrum from a gapless one to one which has a 
finite energy gap as k —> 0. 

As we mentioned above, although it may appear at 
first glance that the wavefunction given in Eq. (|141|) 
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should allow us to reproduce the salient features of Bo- 
goliubov's theory, a careful analysis shows that, in fact, 
different results are obtained for the ground state energy 
and for the energy of elementary excitations, since now 
the number of particles in the condensate is not given 
by (TV — 2m) as was the case in the variational approach 
for the single-mode Hamiltonian , but is given by the 
more complicated expression (N — 2^2iLi n i) involving 
all momentum modes. In particular, in Appendix |B] we 
show that the expectation value of the Hamiltonian 
in the state |\P(iV)) is no longer given by Eq. (|67|) , but 
by the following expression: 



^(N)\H k M(N)) 



where v(kj) is given by: 



£k, 



v(kj)n B 



v{kj)n B - 



C k , 



M 



v(kj) 



k, 



1-ci 



(144) 



(145) 



If it were not for the term between parenthesis in this 
last equation, the result in Eq. (|144[) would be perfectly 
identical to the expectation value obtained within the 
single- mode approach, Eq. (|67l) . Minimization of the ex- 
pectation value of the total Hamiltonian {^(N)\H\^(N)) 
with respect to the variational parameters {ck} would 
then give the same values for these parameters as in the 
single- mode case, Eq. (|7Tj) . and hence everything we dis- 
cussed in Sec. [V] would remain pretty much unchanged. 
Below we will show that minimization of (^f(N)\H\^f(N)) 
as given by Eq. (|144l) . i.e. with ?;(k) replaced by v(k), 
leads to a very different solution for the Ck's, leading 
in turn to some important changes for the ground state 
properties of the interacting Bose system. 

Going back to the result of Eq. (11441) . if we now 
minimize the expectation value of the total Hamiltonian 



H = ^2 j=i > w hich is given by: 



(*(N)\H\*(N)) = 



E 



:{h, 



n B v(kj] 



Ck, 



n B v(kj) \ 



(146) 



with respect to the constants c-^. , instead of Eq. (155)) we 
now obtain the following equation (see Appendix [Cj : 



<- 2 tobr K + 1 = 



(147) 



v n B v(k 3/ 

In the above equation, £kj now denotes the quantity: 

<?kj = £kj + n B v(kj) + cr k , (148a) 



2 M 

with er kj = — n B v{ki)- 

i=i(#j) 



Ck, 



(148b) 



Solving Eq. (| 147|) for Ck, we obtain: 



Ck 



( gk ) 

\n B v(k) I 



<n B v(k) 



- 1. 



(149) 



where again the sign of the second term has been chosen 
so that < Ck < 1. Expressions (|146l) for the ground 
state expectation value of the Hamiltonian and (I149p for 
the coefficients of the ground state wavefunction are the 
main results of this paper. In the rest of this Section we 
want to explore how these expressions, which, to the best 
of the author's knowledge, have not been studied previ- 
ously, alter the description of interacting Bose systems 
given in the standard Bogoliubov formulation. 

We now rewrite the quantity Ok of Eq. (|148b[) in the 
form: 

2 ^ ,, s c kl 2 fu ^ c kj 



l + c kj 

(150) 

Using the fact that n B — N/V, and neglecting the second 
term on the rhs of this last equation, (which is of order 
0(1/ N)) we obtain: 



M 
i=l 



f 'k, 



c k, 



(151) 



Transforming the sum into an integral, we obtain (notice 
that the factor of 2 disappears because the summation in 
Eq. (|15ip is over half of phase space only - notice also 
that the number of momentum modes M has been sent 
to infinity): 



Ok 



Ck' 



C k ' 



(152) 



As it can be seen, the quantity on the rhs of the above 
equation does not depend on k, and we shall henceforth 
drop the subscript k from Ok, and rewrite Eq. (|148b|) in 
the form: 



£kj = £k 7 + n B v(kj) + a, 

with a = I v(k) — - — . 

A 1 + Ck 



(153a) 
(153b) 



In the particular case where v(k) = g 7 corresponding 
to v(r) = gS(r) in real space, it is not difficult to see that 
the integral in Eq. (|153b>[) has an ultraviolet divergence in 
three dimensions, since Ck as given in Eq. (|149l) behaves 
like 1/fc 2 as k — > oo. To circumvent this difficulty, instead 
of the interaction potential v(r) = gS(r) we shall use: 



v(r) 



ge 



-r 2 /(2A 2 ) 



(2ttA 2 ) 3 /2 



(154) 



where A is a positive quantity having the dimensions of 
length. This expression of the interaction potential has 
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been chosen to ensure that we recover the form v(r) = 
gS(r) in the limit A — > 0. In Fourier space, the interaction 
potential of Eq. (|154j) is given by: 



u(k) = ge 



(155) 



With this expression of v(k), the integral on the rhs of 
Eq. (I153b[) converges. Going back to this last equation, 
we see that the integral on the rhs depends on a through 
Ck, and hence we see that Eq. (|153b[) can be seen as a 
self-consistent equation for a. In the following, we shall 
solve this self-consitent equation and find the value of a 
for a given value of the paramters g, ub and A. But, 
before we do that, we want to go back and re-examine 
the quantity v(k) we introduced in Eq. (j!45[) . 
Let us rewrite the rhs of Eq. (|145[) in the form: 



2 A d 



2 

, < li^r)- (156) 



4 N 1 - c k 



The last term between brackets in the above equation 
is of order 0(1/N), and hence can be neglected in the 
thermodynamic N — > oo limit. On the other hand, the 
sum X)j=i c ki/(l — c ki) * s recognized as the expectation 
value (^{N)\Y2^i a L a k i \ i &(N)), where the summation 
extends over half of the wavevectors k in phase space. 
Hence, we can write: 

t;(k)=«(k)(l-^) > (157) 
where iV<j is the total number of depleted bosons: 

N d = 2j2rr i 2- ( 158a ) 

i=l 1 C k, 

= E i^? ( 158b ) 

k#0 1 C k 

For ease of notation, we shall introduce a new symbol Cd, 
which we shall dub the "depletion factor" , such that: 



C d = l 



Nd 
N 



(159) 



Then, Eq. (jl5T[) can be rewritten in the form: 

5(k) = O(k), (160) 

We now are ready to tackle the self-consistent equa- 
tion for a. To this end, we will again use dimen- 
sionless units, where energies are measured in units of 
ubv(0) = griB, and wavevectors are measured in units 
of feo = \/2mnBg/h~i and we will find it convenient to 
express the interaction strength g in terms of the s-wave 
scattering length a, such that: 



4:Trah2 



in 



(161) 



Then we can write for Ck the following expression: 
1 



c k - l + -^(fc 2 + a)e 4 ™ saA fe 
w 



1 

a, 



1 + -^-{k 2 + CT)e 47r "B QA2fc2 



- 1, 



(162) 



where we denote by k and a the dimensionless quantities: 

(163a) 
(163b) 



k 



gn B 



If we use this last expression of Ck in Eq. (|153b[) . and 
change the variable of integration from k to the dimen- 
sionless quantity k, we can write for the dimensionless 
quantity a the following self-consistent equation: 



8\/2, 



dk fc2 e --M™ B aA 2 )fe 2 



1 + Q 2 



Q 2 



where we used the shorthand notation: 



Q 2 = Cj 1 (k 2 + a)e i ^ nBaX2 )~ k2 



(164) 



(165) 



A numerical solution to the above equation for a can be 
found by iteration in the following way. First, one starts 
with an initial guess for a . Using this initial guess, one 
computes the ratio Nd/N using Eq. (|172l) below, which 
allows us to find the depletion factor Cd = 1 — Nd/N. 
This value of Cd is then used to solve the equation (|164|) 
for a . This computed value of a is then used again as an 
input to find a better estimate of the ratio Nd/N using 
Eq. (|172|) . and the process is henceforth repeated until 
convergence and a stable solution for C d and a is found. 

In the following, we shall be mostly interested in dilute 
Bose gases, for which risa 3 <C 1. There indeed seems to 
be a broad consensus in the physics community^ that 
Bogoliubov's theory is not adequate for describing the 
properties of denser systems such as, e.g., superfluid liq- 
uid Helium near zero temperature. For the particular 
case of liquid Helium, this belief stems from the fact 
that Bogoliubov's theory fails to capture a few major 
features of this system, most notably the fact that the 
depletion, even at the lowest temperatures, is quite sub- 
stantive (the number of bosons in the condensate not rep- 
resenting more than 10% of the total number of bosons N 
in the system), as well as the appearance of a roton min- 
imum in the excitation spectrum at higher values of the 
wavevector k. (It can in fact be shown that the failure of 
Bogoliubov's method to describe these specific two fea- 
tures of liquid Helium can be easily remedied through an 
appropriate modification of the Bogoliubov Hamiltonian. 
This discussion, however, is beyond the scope of this pa- 
per, and will be postponed to a future contribution^) 



23 



Given the above, we shall restrict ourselves to the situa- 
tion of a dilute Bose gas, the ground state properties of 
which are believed to be properly captured by the Bo- 
goliubov Hamiltonian. 

For simplicity, we shall take the characteristic length 
scale A which governs the range of the interaction poten- 
tial v(r) to be the scattering length a. Under these cir- 
cumstances, a numerical solution of the self-consistency 
equation (|164[) for n B a 3 — 10~ 3 yields the following val- 
ues for the quantities a and Cd (the reader will note that 
similar results are obtained using other numerical values 
of the parameter nga 3 as well): 



a = 0.3899 
C d = 0.9762 



0.39, 
0.98. 



(166a) 
(166b) 



With the knowledge of a and Cd, we now are in a posi- 
tion to determine the coefficients Ck, and hence find the 
depletion and the ground state energy of the gas. Going 
back to Eq. (|162p . we see that the exponential factor 
exp(fc 2 A 2 /2) for A = a is given by exp(47rnB0 3 fc 2 ), and 
for small values of k this quantity can be approximated 
by unity (since rt^a 3 = 10~ 3 -C 1). The expression of Ck 
can therefore be approximated by: 



Ck 



_o 2 



- 1. 

(167) 



Notice that this expression of Ck reduces to the expression 
we found before in the single- mode model, Eq. (|71|) . in 
the limit a — and Cd = 1- Fig. [TU] shows a plot of Ck as 
a function of the dimensionless wavevector k for a ~ 0.39 
and Cd — 0.98 (lower curve), and for a = and C d = 1 
(upper curve) . It is seen that the limiting value of Ck as 
k — > when a ^ is not equal to unity, but is given by: 



Ck^O = 1 



C7 1 a 



1 - 0.42. (168) 
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FIG. 10: Plot of the coefficient Ck as a function of the dimen- 
sionless wavevector k. The dashed line shows the coefficient 
Ck as given in the single-mode theory, Eq. (|71[) . The solid line, 
on the other hand, shows the coefficient Ck as obtained from 
the multi-mode variational treatment of the present section, 
Eq. (T149)) . with a = 0.39 and C d = 0.9762. 



depletion does not diverge as k — > 0. Here, by contrast, 
Ck goes to a limit which is smaller than 1 as k — > 0, 
and hence iVk does not diverge at small values of k even 
though the geometric summations (over the coefficients 
C n of the wavefunction) in the derivation of Eq. (|169l) are 
extended to infinity. Fig. [TT]shows a plot of the depletion 
iVk as a function of the dimensionless wavevector k for 
the single-mode theory where Ck is given by Eq. (|7ip 
(upper curve), and for the multi-mode approach where 
c k is given by Eq. 1(157]) . with a = 0.39 and C d = 0.98 
(lower curve). 

The total depletion is given by the following expres- 



Nd 



^ i - 

k^O 



(170) 



C. Depletion of the condensate in presence of 
many momentum modes 

With the knowledge of Ck, it is easy to find the number 
of depleted bosons iVk in the single momentum state k, 
which is given by: 



W k = (tf (iV)|4a k |*(iV)) = 



1 



(169) 



In the case of the single mode theory, the coefficient Ck 
as given by Eq. (J7TJ) goes to unity as k — > 0, leading to a 
divergence of the quantity JVk at small wavevectors. This 
happens, as we have shown in detail in Sec. IV Bl when 
we work in the N — > oo limit and extend summations 
over the coefficients of the wavefunction to infinity; but 
otherwise if we insist on working with finite sums, the 



The summand, after a few manipulations, can be written 
in the form: 



l + Q 2 



1 



VOW 2 



(171) 



Transforming the sum in Eq. (|170[) into an integral, 
we can write the ratio Nd/N of the number of depleted 
bosons Nd to the total number of bosons N in the form: 



IF = —f=(n B a )- 



Nd 
N 



dk k z 



1 



(172) 



In the standard formulation of Bogoliubov's theory, 
which corresponds to setting a = 0, Cd = 1 and A = in 
the above expression, the integral in the above equation 
can be evaluated exactly, and has a value of (y/2/3) ~ 
0.4714. For non-zero values of a, it is straightforward 
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FIG. 11 



Plot of the depletion of the condensate iV k = c£/(l- 



c k ) as a function of the dimensionless wavevector k. The 
dashed line shows the depletion obtained using the coefficients 
Ck obtained in the single-mode theory, Eq. (|71[) . The solid 
line, on the other hand, shows the depletion obtained using 
the coefficients Ck obtained from the multi-mode variational 
treatment of the present section, Eq. (|149[) . with a — 0.39 
and C d = 0.9762. 



to evaluate the integral numerically. For example, for 
a = 0.39 and A = 0, the integral evaluates to 0.3276, 
which translates into a reduction in the number of de- 
pleted bosons by about a third (0.3276/0.4714 ~ 0.695) 
with respect to the result of the standard Bogoliubov 
method. 



D. Ground state energy in presence of many 
momentum modes 



We are now in a position to find the ground state en- 
ergy (H) = (^(N)\H\^(N)) of the system, which is given 
by: 



(H) = 



1 ^ c 2 { [ £ k + ^Wfc)]c k - n B v(k)cwY (173) 



Using the expression (I167P of the coefficients Ck into this 
last equation leads to an ultraviolet divergence at large 
k, much like in the standard Bogoliubov approach. To 
circumvent this difficulty, we here shall use the "regu- 
larized" expression (I162[) of Ck, which takes into account 
the fact that the interaction potential between bosons 
falls off at large k. Then, if we transform the sum into 



an integral, we obtain (in three dimensions): 
E 1 



— = --an. 



V 



[b/llnza?) 1 * r dkk 2 

\ V7T Jo 



C' d e 



( 



1 + Q 2 - vWT^ 



1 + Q 2 



(174) 



Numerical evaluation of the above integral for A = a 
and nsa 3 = 0.001 yields, for a = (the corresponding 
value of d, evaluated using Eq. (|172|) . is given by d = 
0.9642), 



^\gn 2 B -(-lM)2)(n B a 3 f^ 

cr—0 Z 

- \gn% ■ (-0.412), 



while for a — 0.39 and Cd — 0.9762, we obtain: 



E 
V 



- lgn 2 B ■ (-U.20)(n B a 3 )^ . 

cr=0.39 z 

\gn% ■ (-0.417). 



(175a) 
(175b) 

(176a) 
(176b) 



We hence see that the solution with a ^ for the coeffi- 
cients Ck leads a lower overall energy than a solution with 
(7 = when the depletion factor Cd is used in the cal- 
culation (we remind the reader that this factor emerges 
when we keep an accurate count of the number of bosons 
in the condensed state k — 0, see Eq. (|145[) ). which con- 
stitutes a direct verification of the validity of our vari- 
ational method. The fact that the solution with a ^ 
has a lower energy than the solution with a = can 
best be visualized in Fig. [T^l where we plot the quantity 
k ("J (N)\Hk\iS> (N)) that appears when we integrate over 
momentum modes in the calculation of the total energy 
E. It appears that the area enclosed by the curve with 
a = 0.39 is greater than the area enclosed by the curve 
with a = 0, hence showing that the value a = 0.39 leads 
to a lower ground state energy of the system. 

One may wonder at this point if the above results for 
the ground state energy are not simply an artifact of 
the special choice we made for the interaction poten- 
tial between bosons, which has a Gaussian dependence 
on k in Fourier space. To clarify this point, in the up- 
per panel of Fig. Q2] below we plot the ground state 
expectation value of as a function of the wavevec- 
tor k from Eq. (|144j) using A = (corresponding to 
v(r) = gS(r)) and with a different value of the parameter 
tibci 3 , nsa 3 = 0.01 (more precisely, we plot the prod- 
uct k 2 (H k ) = fc 2 (*(A^)| J ff k |*(iV)) vs. k, the factor k 2 
representing the Jacobian of the integral over wavevec- 
tors in three dimensions that appears in the calculation 
of the ground state energy E = (^(N)\H\^(N))). In 
this last figure, the solid line shows the result for k 2 (i? k ) 
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FIG. 12: Plot of the quantity k 2 ($(N)\H k \¥(N)) (where 
the factor k 2 comes from the Jacobian of the integral over 
wavevectors in three dimensions) that appears in the calcu- 
lation of the ground state energy E = {^f(N)\H\^(N)). The 
dashed curve corresponds to a = and Cd = 0.9642, while 
the solid curve corresponds to a = 0.39 and Cd = 0.9762. 



one obtains by plugging the coefficients Ck of the single- 
mode theory, Eq. ([71]), into Eq. (fT44|) . The dotted lines, 
on the other hand, show the results for this same quan- 
tity one obtains using the coefficients from Eq. (|167[) 
for three nonzero values of a. It is again seen that the 
expectation values (i?k) that are obtained using the co- 
efficients Ck from Eq. (|167j) with it ^ are consistently 
lower than the one calculated using the single mode co- 
efficients from Eq. (fTTj). which gives further credence 
to our minimization procedure in which the number of 
bosons in the condensate (N — 2 ^ Uj) is kept through- 
out the calculation. Conversely, in the lower panel of 
Fig. (fLlf we plot the quantity fc 2 (V>k|#k|V'k) from Eq. 
(|57| using the Bogoliubov result (|TTj) for Ck (solid line), 
and using the result (jl 6T[) for various nonzero values of 
a . It is seen that for the single-mode model the standard 
result for Ck (which corresponds to setting a = and 
Cd = 1 in Eq. (|167l) ) is always smaller in energy than 
the result obtained by using a non-zero value of a = in 
the expression of Ck from Eq. (|167[) . We therefore con- 
clude that both our minimizations, in Sec. IV Al and in 
the present Section, produce correct functional forms for 
the constants Ck that minimize the ground state energy 
of the system, with Eq. (|71[) representing the correct 
functional form of Ck when we consider the single-mode 
Hamiltonian (or, equivalently, when we use the Bo- 
goliubov prescription ao ~ y/~N and remove the k = 
mode from the Hilbert space, so that all the single-mode 
Hamiltonians are decoupled), and Eq. (|167l) repre- 
senting the functional form that minimizes the ground 
state energy when we do not use Bogoliubov's prescrip- 
tion and keep the correct expression (N — 2J2i n i) of 
the number of bosons in the k = state throughout the 
calculation. 




FIG. 13: Upper panel: Plot of the quantity 
k 2 (^(N)\Hi,\^(N)) using Eq. (fT44)l for A = and 
hbo 3 = 0.01. The solid curve corresponds to a = and 
using the result for the constants c k from Bogoliubov's 
theory, Eq. (171(1 . The dotted curves are obtained by using 
Eq. ([167)) for c k , with b = 0.2, b = 0.4 and d = 0.6 
from top to bottom respectively. Lower panel: Plot of the 
quantity k 2 (^(N)\H k \^(N)} using Eq. (g>7]) for A = and 
riBfl 3 = 0.01. The solid curve corresponds to a = and c k 
from Bogoliubov's theory, Eq. (|71|) . There are three other 
dashed curves which correspond to using Eq. (|167|) for c k , 
with a — 0.2, a = 0.4 and a = 0.6. These three curves are 
very close to one another, and are all situated above the solid 
line corresponding to the result of Bogolibov's theory. 



E. Elementary excitations in presence of many 
momentum modes 



We now want to examine the elementary excitations of 
the Bose system in presence of many momentum modes. 
Let us again define the operators «k and as follows: 

a k = Uka k a|, + Vka a^ w , (177a) 
a k = Uka k a + WkftQa-k- (177b) 
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It is easy to verify that the action of a kl on the ket 
|^(JV)} gives the following result: 



CO oo 



A I 



\ 



n 1 (N+ 1 - 2^: 



1 



M 



m[N + 2-2^', 

i=l 

M 

x |JV + 1 — 2 t^; m — 1, m; 712,712; 



(178) 

Requiring that akJ^iV)) = gives the condition: 



Cni— 1 



AT- 



■2-2E^im 



(179) 



Cni \ JV + 1 

Approximating the square root with unity, we obtain: 

(180a) 



^ki 
Uk a 



C ni -i 
Pki, 



(180b) 



which is the same condition obtained previously for the 
single- mode case, Eq. (19"8"1) . Proceeding in the same way 
as in Sec. IV Dl we can determine the constants Uk and Vk 
by using Eq. (|180[) and the requirement that the excited 
state a k |\&(./V)) be normalized to unity, i.e. that: 



(y(N)\a k al\V{N)} = 1. 



(181) 



The above condition leads to the following result for the 



quantity u k 



U k 



7k, 



(182) 



where now the quantity 7k is given by (compare with Eq. 
(ED): 



7^ 



1 



(a^a ) - 
1 

(aja ) ' 



a k ak) 



(183a) 
(183b) 



In the above equations, we denote by (• • • ) the expec- 
tation value in the ground state |W(AT)), i.e. (•••) = 
($(N)\ ■ ■ ■ |*(A r )). (In Sec. IV D| 7 k was defined in terms 
of expectation values in the normalized ground state |^k) 
of the single-mode Hamiltonian Hk-) Note also that, in 
going from the first to the second line of Eq. (|111|) we 
again assumed that (a k a k ) <C (aja ), i.e. that the deple- 
tion of the condensate into any single-particle state with 



momentum k is small. Using Eqs. (jl 80[) and (|182l) . we 
now can determine the values of the constants £t k and v k . 
These can be written in the form: 



Uk = 7k"k, Vk = 7k«k, 

1 c k 



Uk 



V 7 ! 



Vk 



V 7 ! 



(184a) 
(184b) 



Using the new result (| 149|) for Ck derived in presence of 
many momentum modes, we finally obtain: 



1 (£k + n B v(k) + a 
Ek 



1 , 



2 1 fSk + n B v(k) + a 
Vir — - 1 



(185a) 
(185b) 



2 V Ek 

where now the spectrum Ek is given by the expression: 
Ek = n B v{kWQ 2 (Q 2 + 2). 



(186) 



We now use the general expression (|A12[) of Hk in 
terms of the operators ctk and a k derived in Appendix 
\X[ which we shall rewrite here for definiteness: 



Hk 



, . L -_. K . «k) - 2Bk] (4 Q k + a -k a -k) 

z 7k k 

[S k (w k + t£) - 2A k ] (alal k + a k a_ k ) 

[^k^k - -BkUk^k] (l ak > a kJ + [ a -k- a -kD 

(Bk«k - ^kUk^k) ([al k , ajj + [a_ k! ak])}, (187) 



where the constants ^4k and Bk are given by (the constant 
rjk is defined in Appendix IX]) : 



w(k) D «(k) 



V 



r 



(188) 



Note that the coherence factors itk and in Eq. (|187|) 
are now given by the expressions in Eq. (I185|) . With 
that in mind, we want to choose a value for rjk that will 
make the quantity [£? k (u k +v k ) — 2A k u k u k ] vanish. Using 
this constraint and Eq. (|185j) . we obtain, after a few 
manipulations: 



Ak = Bk(^ + l 

n B v(k) 



(189) 



Using the definitions of A k and i? k given in Eq. (11881) . 
we can find the expression of the quantity r/k, which is 
now given by: 



1 vCk) ( a 

Vk= Tt ^tx 1 + — 
N v(k) \ ek 



(190) 



We now are in a position to calculate the excitation en- 
ergy associated with Hamiltonian Hk, which is the coef- 
ficient of the quadratic terms a k a k and a_ k a_ k in Eq. 
(|187p (it can indeed be shown that the last three terms 
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on the rhs of this last equation do not contribute to the 
excitation energy, see Appendix [A")l . i.e.: 

AE exc (k) = (*(7V)|a kJ ff k aj c |*(7V)) - (^(N)\H k \^(N)}, 

(191a) 

= 7T2 [ A ^ u l + w k) - 2B k K k « k ] . (191b) 
^7 k 



Using the results (| 149)) for the coefficient c k and (|190p for 
the quantity ?y k , along with the definitions (|185[) of the 
constants u k and v k , we hnd, after a few manipulations: 



1 



1 



AE exc (k) = -E k = -n B v(k)y/Q 2 (Q 2 + 2), (192) 

where we remind the reader that: 

Q 2 = (193a) 
n B v(k) 

= C^(k 2 + a)exp(4irn B a\ 2 k 2 ). (193b) 

We again would like to emphasize that AE exc (k) is the 
excitation energy associated with the Hamiltonian iJ k . 
The full Hamiltonian H — X) k ^o contains an iden- 
tical contribution i?_ k , so that the total energy cost 
A£g°* c (k) to bring the system from its ground state 
|*(A0) to the excited state a\.\^(N)} is given by: 



A£*f c (k) = AE exc (k) + AE exc (-k), 



n B v(k)y/Q 2 (Q 2 + 2). 



(194a) 
(194b) 



The above expression of the excitation energy, after a 
few manipulations, can be written in the form: 



A£^(k) _ 1 

" c d 



n B g 



^(k 2 +d)(k 2 + a + 2C d i 



— Ann f>a\ 2 k 2 



(195) 



Again, given the dilute Bose gas parameters we have con- 
sidered, the exponential exp(— fc 2 A 2 /2) is close to unity 
for the most interesting values of k such that k <C A -1 , 
and hence we can write: 



A^ c (k) 1 



n B g 



— \/(k 2 +a)(k 2 +d + 2C d ). (196) 



Notice that the above expression reduces to the usual 
Bogoliubov spectrum for the single-mode theory in the 
limit <T = and C d = 1. In the case where a ^ 0, 
however, a major feature of the spectrum (|196|) is that it 
displays a finite gap as k — > 0, which is given by: 



A£%' c (k = 0) 
n B g 



= </-(- + 2 



(197) 



For a — 0.39, the numerical value of the gap is given by: 
A^° t c (k = 0)~0.98n B5 , (198) 
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FIG. 14: Solid line: Plot of the excitation spectrum A-EgJ* (k) 
of Eq. (|195|) . Dashed line: Hartree-Fock spectrum, Erf = 
£k + ns«(k) = e k + gn s . 



which is comparable to the gap predicted by Hartree- 
Fock theory, which is given by AE^^k = 0) = n B g, 
and to the gap predicted by Girardeau and Arnowitt a 
long time ago^i using a different method than the one we 
employed in the present study. 

At this point, we would like to caution the reader that 
the above result does not by any means imply that the 
actual excitation spectrum of bosons in a real (experi- 
mental) system has to be gapped. It actually only implies 
that Bogoliubov's theory predicts a finite gap in the k — ¥ 
limit when this theory is formulated within a number- 
conserving framework where an accurate count of the 
number of bosons in the condensate is kept throughout 
the calculation. This author would like to make it clear 
that he is not by any means advocating a gapped exci- 
tation spectrum for bosons, but merely presenting what 
Bogoliubov's theory predicts for this spectrum when the 
calculation is performed in a tightly number-conserving 
fashion. Going back to the Girardeau- Arnowitt number- 
conserving theory 2 ^ mentioned above, which also predicts 
the existence of a finite gap in the excitation spectrum, 
an interesting aspect of this theory is a demonstration 
that nonpairing "triplet" contributions to the Hamilto- 
nian are of the right order of magnitude to cancel the 
excitation energy gap. Takano subsequently showed 33 
that such cancellation does indeed occur, and it may be 
possible that a similar scenario takes place for the vari- 
ational model studied in this paper if additional terms 
are included in the Bogoliubov Hamiltonian. While such 
a scenario cannot, a priori, be ruled out, our results are 
still important on their own, because they show that, 
as far as Bogoliubov's Hamiltonian is concerned, the 
conventional excitation spectrum of Bogoliubov's the- 
ory becomes gapped when the calcualtion is performed 
in a number-conserving fashion. More work is certainly 
needed in order to ascertain under what conditions the 
spectrum of excitation of bosons is indeed gapped or not 
in the k — > limit. 
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VIII. DISCUSSION 

Having explored the results one obtains for the main 
physical observables when the full Hamiltonian H = 
Sk^o is diagonalized taking the conservation of bo- 
son number into account, we now want to present a brief 
summary of our results, and discuss the implications of 
these results on the formulation of Bogoliubov's theory 
of an interacting Bose gas. 



A. Summary of our results 

A major goal of this paper has been to clarify the na- 
ture and meaning of the standard formulation of Bogoli- 
ubov's theory. Our investigation has shown that this 
theory seeks to diagonalize each single-mode Hamilto- 
nian Hk independently. Comparison with the results of 
the exact numerical diagonalization of the single-mode 
Hamiltonian Hk have shown that Bogoliubov's theory 
gives an astonishingly accurate description of the ground 
state energy and excitation spectrum of H^. While in a 
number non-conserving framework, where the k = state 
is removed from the Hilbert space used to describe the 
system, one can write the ground state wavefunction as a 
simple product of the ground state wavefunctions of the 
i?k's, and hence diagonalizing the single-mode Hamilto- 
nians independently from one another makes sense, in 
a number-conserving framework diagonalizing the H^s 
independently is not very helpful, since the ground state 
wavefunction of H — X)k^o ^ k cann °t m this case be 
written as a simple product of the ground state wave- 
functions of the Hk's (the reason this is so is because the 
Hilbert spaces spanned by these single-mode Hamiltoni- 
ans have the k = state in common). The above ob- 
servation calls for a more careful diagonalization method 
where the total Hamiltonian H is diagonalized directly. 
This is what we have attempted to do in this paper by 
using the fully number-conserving trial wavefunction of 
Eq. (|141[) . and finding the variational coefficients Ck in 
exactly the same way as in Ref. |29|. It is actually quite 
on purpose that we have presented the calculation of the 
ground state energy in such a level of detail both for the 
single-mode Hamiltonian thereby reproducing all the 
results of the conventional Bogoliubov theory, and for the 
total Hamiltonian H . Comparing the two methods un- 
derscores the stark differences between a theory that does 
enforce the conservation of the total number of bosons N 
and theories which do not, and shows explicitly that the 
results of Bogoliubov's method cannot be obtained when 
N is conserved between all the momentum modes, but 
only in a single-mode approach where N is conserved for 
a single momentum mode. In the following, we shall sum- 
marize a few salient results of our variational treatment, 
and try to discuss how this treatment improves on the 
standard formulation of Bogoliubov's method. 



1. Divergence of the depletion (a^a^) and of the 
anomalous average {a-^a-k) in the k — > limit — As we 
have seen in Sec. MI Al in the standard, number non- 
conserving version of Bogoliubov's theory, both the de- 
pletion (et^etk) and the anomalous average (akOk) diverge 
in the k — > limit. The use of a number-conserving ap- 
proach, as we explicitly have shown in Sections IV Bl and 
IV CI for the single-mode theory and in Sec. I VII Cl in pres- 
ence of many momentum modes, removes these unphys- 
ical divergences and yields results that are always finite 
for finite N. 

2. Use of the Hamiltonian H instead of H — fj,N to 
describe the energy of the system — In addition to re- 
moving the unphysical divergences from the quantities 
(a^ak) and (akOk) in the k — > limit, the variational 
method has the advantage of restoring to the Hamilto- 
nian H its usual meaning as the operator representing 
the total energy of the system. As we saw in Sec. MI Bl 
this is not the case in the conventional formulation of Bo- 
goliubov's theory, where the role of the energy operator 
is played by the combination H — /xiV, and where the use 
of H gives rise to nonsensical results. 

3. Commutation relations of the «k operators — An 
interesting aspect of our investigation is that we have 
clarified the origin of the commutation relations between 
the excitation operators ak and a k . Indeed, the number- 
conserving version of these operators, where ak is given 
by: 

«k = 7k(ukakao + w kal k a ), (199) 
obeys the following commutation relations: 

[a_k, ak] = [a^ k , a kJ = ll u kVk [« k a k - aL k a -k] ■ 

(200a) 

[a k , a k ] a 7k [ala - u k a k a k + vla^_ k a- k ] . (200b) 

As we have seen in Sec. IV Dl imposing the commutation 
relation [ak,a k ] = 1 is nothing more than a convenient 
way to ensure that the excited state a k |\E , (A^)) is nor- 
malized to unity. Given that the rhs of Eq. (|200b[) is an 
operator and not a c-number, it is not possible to satisfy 
the commutation relation [ak, a k ] = 1 for all possible 
states in the Hilbert space, and hence for weak perturba- 
tions one simply requires that this commutation relation 
be satisfied in an averaged sense at the ground state, i.e. 
(*(7V)|[a k ,a k ]|*(iV)) = 1. This led us to the following 
condition on the constant 7k: 



( a a 0) - « a k) 

with the consequence that itk and i>k verify the identity 
w k — v\ = 1, much like in Bogoliubov's approach. It is 
to be noted that, even in the single-mode theory of Sec. 
IV Dl the condition (I20ip by itself is not sufficient to re- 
cover the Bogoliubov spectrum of excitations, for it can 
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be verified that we also require that the depletion of the 
ground state into any given momentum mode k be small, 
(ajao) 3> (a^dk), so that we may write j 2 . ~ l/(aJao). 
Should a situation arise where the depletion of the ground 
state is no longer small (such as for liquid Helium at very 
low temperatures, where the condensate fraction repre- 
sents only about 10% of the total system^), then the full 
expression of 7^ given above, Eq. (|20ip . has to be used, 
and there is then no guarantee that the Bogoliubov spec- 
trum will be recovered, even in the single-mode theory of 
Section |VDj 

4- Single-mode theory: number-conserving variational 
approach for the single-mode Hamiltonian — As men- 
tioned in the opening paragraph to this subsection, in 
this paper we have shown that Bogoliubov's theory cor- 
responds to a decoupled approach in which each single- 
mode Hamiltonian iJk is diagonalized separately from the 
other momentum contributions Hjc'f/k)- What is more, 
we have performed an exact numerical diagonalization 
of the Hamiltonian H k . Comparison of this exact nu- 
merical diagonalization and the variational Bogoliubov 
treatment shows that Bogoliubov's theory gives spectac- 
ularly accurate results for the ground state energy and 
the excitation spectra of each of the single-mode Hamil- 
tonians H^. However, the results of Bogoliubov's method 
for the depletion of the condensate are less accurate, as 
this method, even in its number-conserving incarnation, 
overestimates the depletion of the condensate by about 
one order of magnitude for small values of the wavevector 
k. 

5. Multi-mode theory: variational approach for the full 
Hamiltonian H — The most important result of this pa- 
per, though, has to do with our variational treatment of 
the full Hamiltonian H, instead of the decoupled kind of 
treatment done in Bogoliubov's method where each mo- 
mentum contribution is diagonalized separately. In 
our variational approach, the k = state is restored to 
the Hilbert space used to describe the system (by contrast 
to Ref. 7 where ao is replaced by V~N), and an accurate 
count is kept of the number of bosons in the k = state. 
As a result, we have shown that the coefficients Ck which 
determine the ground state wavefunction of the sys- 
tem are no longer given by the expression (I71[) obtained 
within the standard Bogoliubov approach. Instead, these 
coefficients are now given by the alternate expression 
(|149|) , which has the profound consequence of giving rise 
to a gap in the excitation spectrum of bosons as k — > 0, by 
contrast to the excitation spectrum in the standard Bo- 
goliubov theory which is gapless in that limit. Note that, 
since the geometrical ansatz C n — y/l — c 2 l (~C] i ) n for the 
coefficients C n of the normalized ground state wavefunc- 
tion |V>k) = J2n=oC n \n) gave such accurate results for 
the ground state energy of the single-mode Hamiltonian 
.ffk, we expect the same ansatz to give equally good re- 
sults when the variational method is applied to the total 
Hamiltonian H = X^k^o ^k, as we did in this paper. 
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FIG. 15: Plot of the product k 2 x (V , k|i?k|V>k) vs. k, as ob- 
tained by exact numerical diagonalization of the single-mode 
Hamiltonian for TV = 200 bosons. Here the wavector k 
takes the values k = 0.01, 0.02, . . . , 10. The fact that the prod- 
uct k 2 x (^kl-ffklV'k} in dimensionless units goes to a constant 
—0.25 at large values of k indicates that {ipk\Hk\ipk) goes to 
zero like — 0.25/fc 2 as k — > 00, hence confirming the result of 
the analytical Bogoliubov theory for that quantity. 



B. A comment on the leading correction to the 
Gross-Pitaevskii value of the ground state energy of 
an interacting Bose gas 

Several comments are in order concerning the explicit 
value of the ground state energy of an interacting Bose 
gas. In the standard Bogoliubov approach, it is claimed 
that this quantity is given by the expression (we now 
include the Gross-Pitaevskii term N(N — l)v(Q)/2V ~ 
N 2 v(0)/2V in our equations): 

E Bog = \vn 2 B v{Q) - i J2 ( £ k + "fl«(k) - E k ). (202) 

k/O 

Of course, we now know that the above expression is 
only valid in a number non-conserving approach, and 
is not the correct expression of the ground state en- 
ergy when the conservation of boson number is taken 
into account. However, for the sake of argument, let 
us momentarily ignore this conceptual difficulty, and re- 
view how the final expression of the ground state en- 
ergy for the total Hamiltonian H is calculated from Eq. 
(|202p . In the case where the Fourier components of the 
interaction potential w(k) are all given by a single con- 
stant, v (k) = g, corresponding to an interaction poten- 
tial which is a delta function v(r) = gS(r) in real space, 
the summand in the above expression has the asymptotic 
form —n B g 2 /{Ae\ l ) = —mn 2 3 g 2 /(2h 2 k 2 ), and hence the 
sum diverges in three dimensions. To take care of this 
divergence, a procedure is devised 2 ^ whereby the cou- 
pling constant g is shifted by the infinite (!) quantity 
—mg 2 J (T^fc 2 ) -1 , and to compensate for this shift 
a quantity + / d 3 k.(2n) 3 mn%g 2 / (2h 2 k 2 ) is added to the 
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expression of Eq, with the result: 
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It turns out that the two terms between parentheses on 
the first line of the above equation represent the first two 
terms in the expansion of the s-wave scattering length a 
in terms of g: 
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h 2 J (2tt) 3 k 2 



(204) 



As it stands, the above expression, literally speaking, pre- 
dicts a negatively infinite value of the scattering length a 
(unless, of course, an ultraviolet cut-off is imposed on the 
momentum integral in such a way that the second term 
on the rhs of Eq. (|204[) is smaller than the first, in which 
case a remains finite and positive; note, however, that if 
such a cut-off is imposed, then we are no longer dealing 
with a delta function potential in real space v(r) = gS(r), 
but with an approximate form of the latter). The result 
of the above manipulations is that the integral in the sec- 
ond line of Eq. (I203P is now convergent in three dimen- 
sions, and leads to the celebrated result that is quoted in 
all standard textbooks on many-particle systems, where 
the final result is expressed not in terms of the original in- 
teraction strength <?, but in terms of the scattering length 
a: 
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Strictly speaking, the leading correction on the rhs of the 
above equation does not involve the scattering length a 
defined in Eq. (|204p . but the "bare" counterpart ao such 
that 4iraoh 2 /m = g. However, in the limit of a dilute 
Bose gas where nsa <C 1, one casually replaces ao by 
a (a very questionable replacement, indeed, given that 
the expression of a contains a divergent integral, see Eq. 
dm) ), leading to: 
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At this point, we want to argue that the manipu- 
lations leading to Eq. (|206l) are mathematically un- 
tenable, for what we simply did is shift the divergence 
from the sum of the expectation values X)k^o(^ k ) k 
j X)k^o(^k — £k ^ n B g) to the scattering length a of Eq. 
(|204| . making the latter a divergent quantity, but the 
original divergence of the ground state energy has not 
been actually removed. To clarify what we mean by the 
above statement, let us give the bare interaction strength 
g a numerical value, say g — 1.0 x 10~ 22 J • A 3 (this is 
not a value that is experimentally relevant to any actual 
system, we just made it up for the sake of argument), 



and ask what is the numerical value, in Joules, of the 
ground state energy per particle. The presence of the di- 
vergent integral on the rhs of Eq. (|204j) prevents us from 
giving an answer to this question, and that is what we 
mean by the statement that the divergence of the ground 
state energy has not been removed by the mathematically 
questionable manipulations leading to Eq. (|20G|) . 

In fact, the above manipulations would have been per- 
fectly legitimate had the integral on the rhs of Eq. (|204[) 
been finite and perturbatively small compared to the bare 
interaction strength g. The fact that this is not the case, 
and that the integral we just mentioned is actually di- 
vergent, makes going from ao in the correction term on 
the rhs of Eq. (I205[) to a in Eq. (|206[) very questionable. 
Furthermore, and as we mentioned above, including the 
divergent integral in the definition of the renormalized 
scattering length a, Eq. (|204[) . merely shifts the diver- 
gence to this last quantity. In terms of the bare scattering 
length ao, or the bare interaction strength g, the ground 
state energy is still a divergent quantity (since a itself is 
infinite according to Eq. (1204p ). and hence we see that, 
strictly speaking, the manipulations leading to Eq. (|206j) 
did not actually remove the original divergence of the 
ground state energy of the system from Eq. (|202j) . 

We now would like to pause for a moment, and plot the 

product fc 2 (Hk)k/«s.9 = fc 2 (V'k|-ffk|'0k)/"-,B<7 as obtained 
from the exact numerical diagonalization of the Hamil- 
tonian as a function of the wavevector k. This plot, 
given in Fig. [151 shows that the above product goes to 
—0.25 at large values of k, implying that the exact nu- 
merical result for (-Hk)k goes to zero like — n 2 B g 2 /{Ae\^) 
(in dimensional units) as k —¥ oo, in agreement with the 
result of Bogoliubov's theory. Unlike divergences that ap- 
pear in other physical theories (e.g. the divergence of the 
perturbative expansion of <^> 4 -type models) which can be 
shown to be artificial divergences and therefore need to 
be removed, we here are in the opposite situation where 
the divergence of the sum X)k^o(^ k )k wnen v ( r ) = 9${ r ) 
is not artificial at all, but on the contrary is part and par- 
cel of the exact diagonalization of the Hamiltonians H^. 
Hence, this is a divergence which should not be tampered 
with and cannot be removed because it corresponds to 
the correct behavior of the quantity X)k^o(^k)k when a 
delta function interaction between bosons is assumed. 

With that said, the whole discussion above is some- 
what irrelevant, because the correct ground state en- 
ergy of an interacting Bose gas, when the full Hamil- 
tonian H is diagonalized properly, is not given by Eq. 
(|202[) . Indeed, in the improved variational procedure 
of Sec. IVII1 the ground state energy does not diverge 
like — n B g 2 / (ie^) = — 0.25n 2 B g 2 /e^, but rather like ~ 
—0.25Cf l n 2 B g 2 /e] i (see Fig. [HI), which makes the whole 
procedure described above (consisting of adding and sub- 
stracting the quantity — n 2 B g 2 / (Ae^)) of no great help. 
The root cause of the divergence of the ground state 
energy being the delta function form of the interaction 
potential, one should really only consider more realis- 
tic interactions whose Fourier transform falls off at high 
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FIG. 16: Plot of the product P x (*|ff k |*) vs. k, for a = 
0.39 and Cd = 0.9762 (solid line) when we set A = in the 
interaction potential (|155p between bosons. The fact that the 
product P x (>I'|_&k| > I') goes to a constant that is higher than 
—0.25 at large values of k indicates that (^'j_fik! , I'} does not 
go to zero like — 0.25/P as k 2> 1. 

momentum, hence making the integral giving the ground 
state energy convergent. In this case, the correction to 
the Gross-Pitaevskii energy v(0)N(N — 1)/2V, expressed 
in terms of bare coupling constants, will be negative, as 
it should, and not confusingly 35 positive as in Eq. (|206|) . 

C. Do the operators a k really create collective 
sound waves ? 

We now turn our attention to the a k operators, and 
discuss the validity of interpreting these operators as cre- 
ation operators for collective phonon modes in view of 
the fact that the spectrum of these excitations, when 
the Hamiltonian is properly diagonalized in a number- 
conserving framework, is no longer gapless as is the case 
in the standard formulation of Bogoliubov's method. To 
gain a little more perspective, let us remind ourselves of 
how phonon modes are defined in the case of a periodic 
monoatomic crystal. If we denote by d(x) the second- 
quantized operator describing the displacement of atoms 
from their equilibrium locations, then we can write: 

aw«£(53f)*;[*' fa -<4>-~1. « 

k 

where Wk = ck (c being the speed of sound in the crystal), 
and where a k and ak create and annihilate a phonon 
mode of wavevector k, respectively. We thus see that 
these operators are related to the displacement field of 
the atoms, and it would be rather awkward to try to 
relate the a^s to the operators which describe adding 
or removing one atom from the system. Yet, this is ex- 
actly how the ak's are defined in Bogoliubov's theory, 
where these operators are usually interpreted in terms 
of phonon modes, even though they are not related to 



any displacement field, but merely describe adding or re- 
moving a boson from the system. In the two paragraphs 
that follow, we want to argue that, whether it be from 
the perspective of the standard, number non-conserving 
formulation of Bogoliubov's theory, or from the perspec- 
tive of a number-conserving approach, the a^'s do not 
represent phonons, but correspond to single-particle ex- 
citations instead. 



1. Nature of the a^'s in the standard, number 
non- conserving Bogoliubov approach 

Let us discuss the nature of the excited modes 
and a k that are defined in the standard, number non- 
conserving formulation of Bogoliubov's theory. For defi- 
niteness, let us rewrite the expressions of these operators 
in terms of the boson creation and annihilation operators: 

a k = u k a k + «kaL k , a fc = Wka fc + UkO-k- (208) 

From the above equations, we see that a k annihilates a 
single boson of momentum — k with probability i>k, and 
creates a single boson of momentum k with probability 
itk- As we already mentioned above, given that a k does 
not act on the condensed bosons, nor does it act on the 
depleted bosons with momentum k' ^ k, it is very sur- 
prising that the state o^lv^s) has come to be interpreted 
as a phonon, the latter being a macroscopic excitation 
involving all the bosons in the system. As the expression 
of ak and Qfk clearly shows, these operators represent 
physical processes involving at most two bosons, and can 
therefore hardly be described as "collective excitations" , 
or "phonons" . Another way to see this is to consider the 
following excited wavefunctions: 

|*<T>) - ^H**>, (209a) 
yrv. 

|*L n) >=c£ B) (4) B l*B>- (209b) 

Conventional wisdom teaches us that these are two very 
different states: ) being a state with n "collective 

phonon modes" , while the state of Eq. (|209b[) rep- 

resents a state where n bosons have been added to the 
Bogoliubov ground state. A closer look, however, reveals 
that these two states are actually identical, and corre- 
spond to the same quantum state obtained by adding n 
bosons of momentum k to the Bogoliubov ground state. 
Indeed, using the fact that a k = Mk(*k — ^ka_k, and the 
fact that the operators a k and a_k commute, we can 
write: 

= d n) [u k al - « k a-k]>B>, (210a) 
= c L n) E I ( nl \ i (^k4) m (-^a-k)"- m |* B ). 

m=0 y ' 

(210b) 
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Now, given that | *S?b) = for all values of the wavevec- 
tor k, we see that all terms in the above summation van- 
ish, except for the term m = n, which gives: 

i*L n) > = 4 n) («*4) n i*B>- ( 2U ) 

Requiring that |^/ k ) ^ e normalized to unity results in 

in] 

the following value of the constant C k : 

Ct ] = -L=, (212) 
upon which Eq. (12 1 1[) becomes: 

|*L B> > = ^?I*b> = l4 n) >- (213) 

This proves our earlier claim that the states |^ . ) and 

l^k ) correspond to the same quantum state. Eq. (I213P 
is an important result, which shows that a state of the 

form |$ k l ') — o^I^b) — ^■ a kl^s) ^ oes n °t represent a 
collective phonon mode at all, and represents merely a 
state with one extra boson of momentum k added to the 
Bogoliubov ground state. 

2. Nature of the ak 's in the variational, number- conserving 
Bogoliubov formulation 

We now want to discuss the nature of the au excita- 
tions from the point of view of the variational formulation 
of Bogoliubov's theory discussed in Section lVIIBI In this 
formulation, is defined by Eq. (|177[> . and we therefore 
have for 4 the following expression; 2 ^ 

4 = ""k4 a o + Uka_ k aJ. (214) 

As we already discussed in Sec. IV Dl the two terms on 
the rhs of the above equation represent the two different 
ways the system can be excited from the ground state 
Ivt'(iV)), which is an eigenstate of the total momentum 
operator P with eigenvalue 0, to a state with momentum 
+k. Again here, we can express a k ao in terms of the 
ak's, with the result (see Eq. (|A6[) of the Appendix): 

a t a o = — [ u k4 - Wka-k] • (215) 

Let us apply the operator on the Ihs of the above equation 
to the ground state |\&(iV)) of the operator H . We have: 

4ao|#(A0) = — h<4 - v k a-k] MN)), 

= ^4|*(JV)), (216) 
7k 

where, in going from the first to the second line, we used 
the fact that a-u\^(N)) — 0. As it can be seen from 



this last equation, the quantum state aj c |^E'(iV)) is iden- 
tical to the state aj c ao| v E'(A)). Since the latter represents 
a state where one boson has been removed from the con- 
densate and a boson has been added to the single-particle 
state with momentum k, we conclude here again that the 
operator aj., when applied to the ground state |^(iV)), 
does not create a collective sound wave at all (such as it 
would be created by shaking the walls of the container, 
or varying the parameters of the confining potential in 
a trapped system — that is after all how a sound wave 
would be created experimentally^ — in which case all 
modes of the system would be excited in a collective way, 
not just one momentum mode as is the case here), and 
merely promotes a single boson from the condensate with 
single-particle momentum k = to the state with single- 
particle momentum +k. 

To summarize this subsection, we conclude that the 
gapless phonon-like spectrum of the 4 operator in the 
standard Bogoliubov method, which in effect creates ex- 
cited states for the single- mode Hamiltonian i? k , has no 
particular physical meaning in terms of travelling sound 
waves. Successive application of the operator 4 to the 
ground state |"0 k ) of the Hamiltonian merely trans- 
fers bosons from the condensate to the single-particle 
state e* k ' r /VV with momentum +k, and hence the cor- 
responding excitations are in no way collective, and can- 
not be identified as sound waves, as they routinely are 
in the literature. This conclusion is corroborated by the 
results we obtained in Sec. IVII El where we have shown 
that, when the k = state is restored into the Hilbert 
space used to describe the system and the conservation of 
the number of bosons is properly taken into account, the 
spectrum of the a k operators has a gap in the long wave- 
length limit k — > and cannot therefore be interpreted 
as a propagating density disturbance corresponding to a 
collective sound wave. (Incidentally, it is quite instruc- 
tive to recall that in the BCS theory of superconductivity, 
the elementary excitation operators are defined in a way 
that is very similar to how the a^'s are defined in Bo- 
goliubov's method. Yet, in BCS theory, the a k are not 
described as phonons, but are instead correctly identified 
as single-particle excitations.) 

D. Comment on the apparent violation of 
Hugenholtz-Pines theorem 

We now want to address the criticism that will in- 
evitably be made about the fact that the variational the- 
ory presented in this paper predicts an excitation spec- 
trum which seems to violate the celebrated Hugenholtz- 
Pines theorem* (HPT). In its simplest expression, this 
theorem states that the zero wavevector and zero fre- 
quency limit of the diagonal and off-diagonal self- 
energies, Sn(k, uS) and £12 (k, lj) respectively, are related 
to the chemical potential fx through the equation: 36 

fiSii(0,0)-fiSi 2 (0,0)=/i. (217) 
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In the standard field-theoretic formulation of Bogoli- 
ubov's theory an immediate consequence of this theorem 
is the emergence of an excitation spectrum which does 
not have a gap as k — > 0. In the formalism presented in 
this paper, however, a finite gap is found, which can be 
seen as an indication that our variational theory violates 
the HPT. 

As an answer to this potential criticism, we shall note 
that the HPT has only been established in the num- 
ber non-conserving field-theoretic formulation of Bogoli- 
ubov's theory. Whether this theorem holds in number- 
conserving situations is far from obvious, and is not, by 
any means, guaranteed^ It is in fact perfectly conceiv- 
able that a number-conserving field-theoretic formulation 
of the theory of interacting bosons may not satisfy the 
same Hugenholtz-Pines theorem satisfied by the number 
non-conserving version, in such a way that the excitation 
spectrum resulting from the modified Hugenholtz-Pines 
theorem for the number-conserving theory is not inconsis- 
tent with the excitation spectrum derived in the present 
study. To construct a number-conserving field theory 
for interacting bosons is, however, a nontrivial exercise. 
What is more, to obtain an equivalent to the Hugenholtz- 
Pines theorem in such a theory may not even be possible 
in the first place, given that a number-conserving the- 
ory would most likely not require the introduction of the 
off-diagonal self energy £12 (k, uS) which appears in the 
HPT. Making any further speculation about this issue is 
far beyond the scope of this paper, and we shall therefore 
content ourselves with the observation made at the begin- 
ning of this paragraph about the fact that the HPT was 
established within a number non-conserving framework. 
As a result, any claim about our variational theory be- 
ing in violation of the HPT would be totally unfounded, 
since this theory is a number-conserving one, and we sim- 
ply have no idea how one can write down an appropriate 
generalization of the HPT that would be valid in number- 
conserving situations. 



E. Consequence on field-theoretic formulations of 
Bogoliubov's theory 

We now would like to make a brief comment on the 
field-theoretic formulations of Bogoliubov's method, both 
at T = and at finite temperatures. It is important to re- 
alize that the field-theoretic formulations currently in use 
in the literature (see for example Ref. [l^ and references 
therein) are based on the Bogoliubov ground state, which 
is a simple product of the individual ground states of the 
single- mode Hamiltonians H^, as in Eq. (|134l) . Needless 
to say, the Green's and spectral density functions will 
have a completely different structure when the field the- 
ory is formulated using the ground state discussed in this 
paper as a starting point. As a result, various other quan- 
tities, such as density-density correlation functions^ and 
the damping of quasiparticle excitations^ may end up 
having expressions that are qualitatively very different 



from those derived within the standard, number non- 
conserving Bogoliubov approximation. 

With that said, it would also be interesting to probe 
whether our variational ground state can be reproduced 
using path integration methods. Path integrals have 
emerged as an elegant and powerful tool to study many- 
body systems, and a successful formulation in terms of 
path integrals which is able to describe the ground state 
studied in this paper may help put the Green's function 
formalism of interacting bosons at finite temperatures on 
a firmer ground. 38 However, such a path-integral formu- 
lation may technically prove to be difficult to achieve, 
since it would most likely require going beyond standard 
Gaussian integration. 

F. Recent literature on number-conserving 
formulations of Bogoliubov's theory 

Before we conclude, we want to briefly discuss recent 
attempts to overcome the difficulties that arise due to 
the non-conservation of particle number in Bogoliubov's 
theory, as our paper would be rather incomplete without 
such a discussion. One notable such attempt is the in- 
teresting work of Castin and Dum, Ref. |25|, where the 
ground state of the interacting Bose system is found by 
writing the field operator ip(r) as a sum of a condensate 
contribution 0o(r)ao and a contribution 8^> (r) from the 
k ^ states: 

^(r) = </» (r)a + ^(r), (218) 

where 0o( r ) denotes the wavefunction of the condensed 
bosons. Then, the quantity <fa/>(r) is treated as a pertur- 
bation with respect to the condensate part </>o(r)ao, an 
approximation which is expected to be valid for a dilute 
Bose gas. Imposing a condition of the form: 

ak|*)=Q, (219) 

which is reminiscent of Eq. (|16[) , and performing a sys- 
tematic expansion in the small parameter e = W Nd /No 
(Nd being the total number of depleted bosons) allows 
one to recover the results of Bogoliubov's theory for the 
depletion of the condensate and the excitation spectrum 
of the system. Overall, the approach presented in this 
work is quite thoughtful, and has the advantage of provid- 
ing a clean derivation of the standard Bogoliubov results 
that overcomes the need to break the U(l) symmetry, 
hence showing that this universally accepted paradigm is 
not by any means required to describe condensed Bose 
systems. By trying to keep the conservation of the to- 
tal number of bosons intact, it avoids many pitfalls of the 
conventional Bogoliubov method. There is a reason, how- 
ever, why this approach yields Bogoliubov-type results 
instead of the results found in the present study. This 
reason has to do with a number of approximations that 
are made at a few key steps of the calculation, and more 
specifically the approximation which consists in replacing 
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the combination a ao by the total number of bosons N. 
This approximation being, to a certain extent, equivalent 
to Bogoliubov's prescription, it is not at all surprising 
that Bogoliubov's results are recovered. By contrast, in 
our variational approach, no such approximation is used, 
and a great deal of care is exercised in order to keep track 
of the exact number of bosons (N — 2 ^T, n,) in the con- 
densate (see Appendix[Bj . The sharp distinction between 
using the approximation ajao ~ N on one hand, and us- 
ing ajao = (N—2 n,i) on the other hand, is what leads 
to the differing results between the two methods. 39 

We now want to comment on the variational approach 
used by Leggett in Ref. |29j, an approach which has been a 
major inspiration for the present study. In this approach, 
the ground state wavefunction of the system is given by 
(see Eq. (8.9) of Ref. M): 

*at = Nl-^falal - ^alal^ N/2 \0). (220) 

k>0 

Although the way this wavefunction is written is dif- 
ferent from the way we wrote the wavefunction in Eq. 
(I14ip . a closer look reveals that these two wavefunc- 
tions are in fact similar to each other in that both rep- 
resent an expansion in terms of the basis states \N — 
2 n i) n ii n i'i ■ ■ ■ j n ii THl ■ ■ •}• I n the course of the calcu- 
lation, however, and in order to find the expectation value 
of the single-mode Hamiltonian in the Bogoliubov 
ground state, Leggett only retains the following part of 
*at (see Eq. (8.15) of Ref. H): 

^ = Nr^(l-ciy/ 2 (alal~c k ataU) |0>. (221) 

This is exactly the ground state wavefunction (l50l) of the 
single-mode Hamiltonian H^, which we call \ipk) in our 
manuscript. By contrast to the wavefunction given in 
Eq. (|220l) above, the variational wavefunction in Eq. 
(I221[) does not take into account the fact that the num- 
ber of bosons in the condensate is given by (N — 2 £V n, ) . 
Hence, this variational calculation uses the same kind of 
approximation we discuss in Sec. |VJ in which the varia- 
tional constants Ck correspond to minimizing the energy 
of each single-mode Hamiltonian independently from 
one another. Again, since the k = state is shared by 
all the wavefunctions \ipk}, the above procedure to find 
the Ck's is mathematically inaccurate, even for a dilute 
Bose gas. This is corroborated by the fact that, when 
the variational calculation is done in a more careful way, 
where the number (N — 2 ^ rii) of bosons in the k = 
is kept throughout the calculation, the results obtained 
for the ground state energy and the excitation spectrum 
are quite different from those of the standard Bogoliubov 
method, as we have explicitly shown in Sec. IVIII of the 
present paper. 

Another attempt at formulating a number-conserving 
theory for the ground state of interacting bosons was 
given by Gardiner in Ref. [23|. We will not discuss this 



study here, but rather point to a comment by Girardeau, 
Ref. [24], where Gardiner's theory is discussed in great 
detail, and is shown to be a special case of the theory 
developped by Girardeau and Arnowitt in 1959, Ref. [2l|. 



IX. CONCLUSIONS 

To summarize, in this paper, we have given a detailed 
and rather thorough discussion of Bogoliubov's theory of 
an interacting Bose gas. Our main point was to reassert a 
result which in principle should be quite well-known but 
is often overlooked in the literature, having to do with 
the fact that Bogoliubov's theory is a theory in which 
each of the sinlge-momentum contributions to the 
total Hamiltonian H = X^k^o ^ k * s diagonalized inde- 
pendently from the other contributions (^k) , and the 
ground state wavefunction of the total Hamiltonian H is 
written as a simple product of the ground state wavefunc- 
tions of the i?k's. As a way to illustrate this point, we 
have discussed a variational, number conserving formu- 
lation of Bogoliubov's method, where the ground state of 
each single-mode Hamiltonian iJk is found independently 
from the ground states of the other Hamiltonians -ffk'^k, 
and we have derived most of the results of Bogoliubov's 
standard, number non-conserving formulation within this 
variational method, including the gapless excitation spec- 
trum predicted by Bogoliubov's method. Arguing that 
the above decoupled way of finding the ground state of 
the total Hamiltonian H may not be accurate, we gen- 
eralized the above mentioned variational method to the 
total Hamiltonian H, making sure to keep an accurate 
count of the number of bosons in the k = state, which 
led us to results that are quantitatively different from 
the results of Bogoliubov's method for the depletion of 
the condensate and the ground state energy of the sys- 
tem. More importantly, it led to an excitation spectrum 
of bosons which has a finite gap as k — > 0, by contrast 
to Bogoliubov's method where this gap vanishes. Given 
the above, the author hopes he has made a compelling 
case that the standard formulation of Bogoliubov's the- 
ory, where not much importance is attached to the con- 
servation of boson number, is far from accurate and hence 
highly unsatisfactory, to say the least, and is in need of 
a major revision. It is the author's hope that the dis- 
cussion presented in this work will help bring about such 
a revision, leading eventually to a more satisfying de- 
scription of dilute Bose systems, both at T — and at 
finite temperatures. For denser systems, such as liquid 
Helium at low temperatures, the Bogoliubov Hamilto- 
nian no longer provides an accurate description of the 
interaction between bosons. A detailed discussion of the 
modifications that need to be brought to Bogoliubov's 
theory to allow us to account for the relatively large de- 
pletion observed experimentally in this system, as well 
as of the roton minimum in the excitation spectrum, will 
be published elsewhere* 3 -^ 
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Appendix A: Diagonalization of the single mode 
Hamiltonian Hk 

In this Appendix, we show how one can diagonalize the 
Hamiltonian Hk in the number-conserving Bogoliubov 
approach. We shall here mostly focus on the case of the 
single-mode theory of Sec. |Vl the general case being very 
similar, apart from a few minor differences discussed at 
the end of this Appendix. 



1. Derivation of the diagonal form of Hk in terms 
of the Qk operators 



Let us rewrite the expression of Hk'. 



i i t t \ v s ) i t 

Hk = ^kla^ak + a^ k a-kj + (a^aoa^ak 

+ t4a aL k a-k + a k a -k a o a o + ao a o a k a -k) ■ (Al) 



We now introduce the operator 6k such that: 

6 k = a,k%. (A2) 

It then follows that the operator «k = M^k^o + ^k a -k a o 
can be written in the form: 

a k = u k 6 k + WkfrLk- (A3) 
Taking the adjoint of the above equation, we obtain: 

ai = u^bl + v k 6_ k . 



(A4) 



These last two expressions of ak and a k in terms of the 
6k 's can easily be inverted to give the expressions of 6k 
in terms of the a k 's, with the result: 



bk = -o [uk^k - vuaV\ . 
7k 



(A5) 



where we remind the reader that 7 k = it k — vl- Using 
the fact that u k = 7k ""k and Vk = 7k Vk, we finally can 
write: 



bk = — [u k a k - «k« k ] 



(A6) 



Now, it is not difficult to see that the pairing terms can 
be easily written in terms of the 6k 's: 

a k; a -k a o a o = a k a o a -k a o = b k b -k' (A7a) 
ajaja k a_k = aja k aja_ k = ^k^-k- (A7b) 

On the other hand, we can write for the Fock term: 



a k Okaga = a k a k (aoaQ - 1), 
= a^aQakO^ — a^ak, 
= b l b * ~ a W- (A8) 



Hence, we can rewrite for Hk the following expression: 

1 ( w ( k )\ r t t \ 

Hk = - [Sk - — J «a k + a!_ k a_ k ) 



2V 



(6 k 6 k + 6 f _ k 6_ k + 6 k 6 f _ k + 6 k 6_ k ) . (A9) 



To the above Hamiltonian, we add and substract the 
quantity hr]k£k{b k bk + 6^ k 6_k)) and rewrite the result 



in the form: 



Hk = H X k + ~(ekV* + ^y-) (6 k 6 k + 6 f _ k 6_ k ) 

+ ^(^tk + 6k6-k), (A10) 
where the excess Hamiltonian H±k is given by: 
Hik = 2 £ k(a k fl k + a^ k a_ k ) 
- ^£kVk(bkbk + &t k &-k) 

-^(4«k + aL k a_ k ). (All) 

For weak perturbation, this Hamiltonian will have a very 
small expectation value, owing to the fact that the two 
terms on the first and second lines cancel each other (we 
remind the reader that 6 k = a k aj, and hence 6 k 6 k = 
a k a k aJao; on the other hand we will see below that rjk = 

7 k , and since 7 k ~ l/No, we see that 7 k 6 k 6 k ~ a k a k ), 
and will henceforth be neglected altogether. (Note that, 
because it only involves a single momentum mode, the 
term on the last line on Eq. (|A11I) can be neglected in 
the thermodynamic limit V — > oo). 

We now can use Eq. (|A6p to rewrite the part of i? k 
that is quadratic in the 6 k 's in terms of the a k 's. A 
somewhat tedious but straightforward calculation gives: 

Hk = -^j [Ak(ul + vl) - 2Bk] (a[a k + al k a_k) 
z lk L 

+ [Bk(ul + vl) - 2Ak] (4«lk + "k«-k) 

+ [Akvl - BkUkVk] ([ak,4l + [ a -k,aL k ]) 

+ {Bkvl - AkUkVk) ([a f _ k , ajj + [a_ k , ajj) }, (A12) 

where we defined the constants A k and Bk such that: 



«(k) „ «(k) 
Ak = ekVk + ^r-, B k - 



V 



V 



(A13) 



We now want to fix the value of the constant rjk in such 
a way as to make the coefficient of the term (a^a^ + 
akCK-k) vanish, i.e. by imposing the condition: 



[Bk(ul + vl) - 2Ak] = 0. 



(A14) 



Using Eqs. (|A13|) and the expressions ([TT2> ([TT31I of the 

coherence factors u k and £) k , we obtain, after a few steps: 



Vk = 7 k - 



(A15) 
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Using the definition of a k in terms of the a k operators, 
one can easily evaluate the commutator [a_ k ,a k ], with 
the result: 

[a_ k , Ok] = [at k , al\ = 7 k WkVk [a k a k - al k a_ k ] . 

(A16) 

For all practical purposes, the last term on the rhs of the 
above equation can be neglected, since for weak pertur- 
bation 7 k ~ 1/iVo ~ 1/N and the expectation values of 
a k <2k and a_ k a_ k should be almost equal, hence cancel- 
ing each other. We therefore can approximate: 



[a_ k ,a k ] = [aL k ,a k ] - 0. 



(A17) 



On the other hand, we have already established in Eq. 
(|108|) of the text that the commutator [afc, ajj is given 
by: 



7k a 0°0 - Wk a lc a k + ^k a -k a -k- ( A18 ) 



Given that S k = 7k"k = Uk/VNo, and v k = JkVk 
v^/y/No, assuming weak depletion the second and last 
terms on the rhs of the above equation can be neglected. 
This implies, since 7 k ~ 1/Nq that we can approximate: 



[a k ,a k ] « 7 k aja ~ 1. 



(A19) 



Under these circumstances, if we use the approximation 
7 k ~ 1/JVb ~ l/N, so that U7 k ~ l/n B , and then use 
the expressions of u k and i> k in terms of c k , Eqs. (|112[) . 
one can show that: 

n B v(k)(ul + vQ = 2(e k + n B v(k))u k v k , (A20) 
then the expression of iJ k reduces to: 

H k = | [e k + n B v(k)] v k - n B v(k)ukVk\ 

^{ [ £ k + n B v(k)] (u k + vQ - 2n s w(k)M k t; k | 



+ 



x (a k a k + aL k «-k)- 



(A21) 



If we again use the expression of u k and v k in terms of e k , 
n B and v(k), Eqs. (|113[) to evaluate the terms between 
curly braces, then we obtain after a few manipulations: 

Hk — ~ [ek + n B v(k) - Ek] + ^£? k (a k O! k + aL k O!- k ), 

(A22) 

which is nothing but Eq. (|116l) of the text. 



2. Proof that the excess Hamiltonian Hik gives a 
negligible contribution to the excitation energy 

In this Subsection, we want to make sure that the ex- 
cess Hamiltonian H\k does not contribute to the excita- 



tion spectrum. We are interested in evaluating the fol- 
lowing quantity: 

Affile = (*(JV)|okffik4|*(JV0) - (*(iV)|#ik|*(iV)), 

(A23a) 

= (^(A0|[a k) #ik]4l*(A0>: ( A2 3b) 

where, in going from the first to the second line, use has 
been made of the commutation relation [a k ,a k ] ~ 1. In 
what follows, we shall find it convenient to decompose 
H\k into the follwing sum: 



with: 



Uk 



H\k — hik + /l2k, 



v(k)\ , * i \ 

k y- J (,a k a k + a7 k a_ k J 



h 2 k = -2 £ kVk( b k b k + b -k b -k) 



(A24) 

(A25a) 
(A25b) 



_ Calculating the commutator [ak,^ik], we find: 



K, hik] = 7X5 ( £ k - — fT - ) [(«k + ^k) a ^ _ 2u k w k aL k ] • 

^k * 

(A26) 



We therefore can write: 



(*(iV)|[a k ,ft lk ]4|*(JV)) = - T if k 

z 7 k 

x (*(N)\(ul+vl)akat\*(N)), 



t>(k) 



w£ + V? 



k -rv k / «(k) 
£k — 



(A27a) 
(A27b) 



27^ v 

Similarly, calculating the commutator [a k , /12k], we find: 

[ofe, h 2 k] = -TrjZkVk [(«k + ^kW - 2tt k v k a^ k ] > 
z 7 k 

(A28) 

and hence: 

(*(^)|[ak^2k]4l v I / W> = -£k??k^4^. (A29) 

7k 



Adding the two contributions from Eqs. (IA27[) and (|A29[) 
together, we find: 



ul + vl 



(*(iV)|[ak,ffik]4l v l / (^)> - " kT ' k 
«(k) 



27^ 



£k 



7k 7 



V 



(A30) 



Given that ?y k = 7 k5 and using the fact that u k = 7kUk 
and «k = 7k «k, we finally obtain: 



(*(JV)|[a k>J ffi k ]al|*(AO) 



■^(4 + 4), (A31) 
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which indeed gives a negligible contribution to the exci- 
tation energy in the thermodynamic limit V — > oo. 

At the end of this Appendix, we briefly discuss the di- 
agonalization of -ff k when a multi-mode model governed 
by the full Hamiltonian H = X^k/o ^ k * s considered. In 
this case, the coherence factors are given by Eq. (|185|) . 
and hence the requirement (|A14|) above gives a different 
result for rj^ than Eq. (|A15|1 . namely Eq. (|190l) . This, 
in turn, leads to the boson excitation spectrum given in 
Eq. (TIM)) of the text. 



Appendix B: Expectation value of the Hamiltonian 
when we keep an accurate tally of the number of 
condensed bosons 

In this Appendix, we give the salient features of how 
we calculate the expectation value of the Hamiltonian 



and hence: 



+ aja a k .ak 3 + a a al. aJ_ k . + ajaja ki a-kj) (Bl) 



in the variational ground state |\1/(A^)) of Eq. (|141l) . 
namely: 



oo oo 



i=l 



\*(N))=ZJ2--- E C m---C nM 

ni— tlm—O 
M 

X |JV — 2^n 4 ;ni,ni; . . .;n M ,n M )- (B2) 

We have: 

<i.a kj MN)) = Z ■■■ E C ni ---C nh 

ni—Q um—0 

M 

x \N -2^n i ;n 1 ,n 1 ;...;nM,n M ), (B3) 



OO OO 



(y(N)\ala k] \V(N)) 



oo 7^ 
n j= ^rij 

r3 



1 



(B4a) 
(B4b) 



where, in going from the first to the second line, we made 
use of the fact that C nj = (—Ck j ) nj , and of the summa- 
tion formulae (valid for |a:| < 1): 



oo 

5>" = T^. 



n=0 
oo 



E-" = r 



(B5a) 



(B5b) 



n=0 



An identical result is obtained for the expectation value 
of the operator al k a_ kj .. 

Let us now turn our attention to the Fock terms. We 
have: 



oo oo 



A I 



m— "Hm— 
M 



x n 



,(Ar-2E^)l^-2E 

i=l i=l 

(B6) 



We shall rewrite the above result in the form: 



oo oo 



M 



^a al.a kj \V(N)) = Z E '" E C ni ■ ■ ■ C nM nj(N - 2nj)\N - 2^n i ;n 1 ,n 1 ; . . .;n M ,n M ) 



m— riM— 

oo oo 



A I 



A I 



■2z E ••• E r " ■■■ ( '» ".•( E ^)i^- 2 E 

ni\nx,nx\ ■ ■ ■ ;n M ,n M )- (B7) 

i=l(jij) i=l 
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Hence: 



oo oo 
X • • • X 
7ii—0 n.j — 

oo oo 



(*(A0l4ao4.a kj MN)) = Z 2 £ x • • • x £ (iV - 2n,) x • • • x £ C 



2 

JIM 



n M =0 



E (E c «r' x E< x - x E< x - x E c »M). ( B8s 



n i= 



nj— 



n M =0 



Z^n,-=0 ^nj 



-2 E 



n»=0 ^Wj A Z^n i =0 < -'n,- 



(B8b) 



Using the summation formulae of Eqs. (|B5[) . and the additional formula 



E 

n=0 



n 2 a; n 



2a; 2 



(l-x) 2 (l-x) s 



we obtain, after a few manipulations: 



(B9) 



(y(N)\alaoal a kj \y(N)) ~ ^t^j- 



\ " 



TV ^ 1 

i=l(#j) 



(BIO) 



Note that, in arriving to this result, a term of the form —2 



has been negelected (this is legitimate 



1-C k 2 (i-c k 2 )3 

given the fact that the terms kept have an overall factor of AT at the front, which makes the term neglected very small 
by comparison in the thermodynamic N — > oo limit). Note also that an identical result is obtained for the expectation 
value of the other Fock term, a aoa_ k .a— k 3 .. 

We now want to find the expectation value of the pairing term aoaoa^ a^ .. We have: 



oo oo 



a a at/_ k .\y(N)} = Z y jr--- ]T C ni ■ ■ ■ C nu ^ {n + If (n - 2^ m) (n - 1 - 2^: 

ni— nM— ) i— 1 i— 1 

M 

x ]jV - 2 - 2 71,; m, m; nj + 1, n 3 - + 1; . ■ .;n M ,n M ), 

i=l 

X. X JW 

£ C ni ---C„ M (n J +l)(jV-2n J ---2 £ 



m— nM— o 



,1/ 



x |JV - 2(nj + 1) - 2 nj;m,ni; . . . ;rc,- + + 1; . . . \n M ,n M ), 

where, in going from the first to the second line, we used the approximation in Eq. ([5§f to write: 



(Blla) 



(Bllb) 



M M - M 

(N-2j2n i )(N-l-2 1 £n i )^N-2n j ---2 £ 



i=l(#j) 



(B12) 
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If we now change the index of summation for the momentum kj from rij to nij = rij + 1, we can write: 

3 



oo oo 



a a 



ni— nij—1 um—0 
A I 

x\N-2rrij-2 ^2 ni;nx,nx;...;mj,mj;...;nM,nM) 

M / oo oo oo 

— I ■•• •■■ C ni C n2 ■ ■ ■ C mj _i • • • C nM mjUi 

i=l(=£j) \ni=0 m 3 =l nj\/=0 
M 

x|iV-2mj-2 2J n i ]n 1 ,n 1 ]...]m j ,m j ;...;nM,n M )- 

i=l(/j) 



(B13) 



Hence, we obtain for the expectation value: 



AT - 2mj + | 



^ E 



Eoo p 2 
rrij=0 rrij 



(B14) 



It can be verified that the only term that one has to keep in the first summation is the term proportional to N (note 
that this is also how the calculation is conducted in the single- mode variational theory of Sec. [V}. Then, if we use 
the summation formulas (|B5p , we obtain, after a few manipulations: 



(*(N)\aoaoal.ai u .\*(N)) = -N 



f'k, 



M 



C k 2 



(B15) 



A similar contribution is also obtained for the expectation value of the operator ao a o a ik a -k 
Collecting all terms, Eqs. (|B4b|) . (IbTU)) and (|B15]) , we obtain: 



(*(JV)|^ Wl |*(JV)) = ek ,-J^-+n fl «(k J )-Jy 1-- £ 7 

-L Oi, . ± Oi, . L iV 1 C 



A/ 



' k 



<n,Bv(\tj) 



c kj 



2 ^ c k2 



? r 



JV — 1 - c: 

i=i(?y) k< 



(B16) 



For the total Hamiltonian H — > we obtain: 



M fr / 9 M c 2 

(*(JV)|^|*(JV))=X;U^+nfl«(ki)(l-^ E 7^7 



2 ^ c, 2 



n fl «(kj)(l-^ E I 



i=l(#) 



AT ^ 1 - c? / 1 - c? 

i=l(^j) k . k i 



(B17) 



which is nothing but Eq. (I146p of the text. 



Appendix C: Minimization of the expectation value over the coefficients c kj . To this end, we shall rewrite 
of the total Hamiltonian this expectation value in the form: 



In this Appendix, we show how we minimize the ex- 
pectation value of the Hamiltonian given in Eq. (|I46[) 
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{9(N)\H\9(N)) 



£k, + n B i 



2 «^ ci 



2 ^ c kj \ c kj 



r 2 y i _ 2 



M fr 2 M 4 



Taking the partial derivative with respect to c k , we obtain: 



2 n B v[k e )(l 

t Ci_ , 



>(*<)(l 



2 A 



iV 



5 1 ) 1 



(CI) 



i=i(/<) 



g(g) 
9c k , 



/ 2 M 4 

£k . +n B u(k i )(l-- ^ — i 
i=l(#j) 
2 2c k , 



2c kj 



«j\ k; 11C,,. 



i=l(#j) 



M 

y3 {"-Wk<? 



*=i(#) 



N 



l-et 



Rearranging the terms on the rhs of the above equation, and using the notation 



2 A c2 



=i(/i) k * 



1 + 



(C2) 



(C3) 



we can write: 



d(H) 
dc kl 



£k, + n B v(kj) -jjY n Bv(ke)—j 



c kf (ck f -l)l 2c kj . 



n B w(kj)- 



1 + 4, 



r 



(C4) 



We hence see that the equation (d(H)/dc kj ) = can Equation (|C5|) can in turn be rewritten in the form: 
be written in the form: 



1 



(1 - 4, 



2£ kj c kj -n B v(k J )(l + cl j ) =0, (C5) 



with 



M 



= Ekj + n B v(k 3 ) + — 53 ™WMy 



c k „ 



£k, 



n s u(kj-) 



c kj +1 = 0, 



(C7) 



(C6) 

which is nothing but Eq. (|147l) of the text. 
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